I don't understand enough about the geometry to visualise the distortions and phase changes that occur off the perpendicular axis. I can see you deal with this in the matlab code. Is that a re-implementation of Ribner and Roys equations? It's hard enough to imagine in a 2D plane, nevermind 3D
All I did with Pd is make a recirculating buffer that has a filter to approximate propagation losses and add N-waves created by the vline segment randomly to it. Sometimes the superposition does create thunder like effects, but that's more luck than judgement it seems.
I reckon to get the right effect you either need several n-wave generators in parallel, or several parallel delays operating on one wave.
As the text says, it's equivalent to a convolution of the n-wave with a set of points that are distances to corners in the tortuous line, it's the hypotenuse of a right angle triange if we assume the bolt comes straight down that's the square root of the horizontal distance squared plus the height squared, and with c as a constant then time correlates directly to distance.
In a way it's granular synthesis, so the density must be calculated. I think you could reduce the whole caper to two variables, density and off-axis phase shift. It's the propagation time divided by the period of one n-wave I think. At ground level the observer is perpendicular so it's very dense, but as you move up the lightning bolt the subtended angle increases and so does the propagation time, so density tails off.
I'm really looking forward to hearing where this takes you next...