Finite difference documentation typo?

[mikessut]

I'm trying to understand finite different propagation and am looking at the documentation here:
https://opticspy.github.io/lightpipes/manual.html#free-space-propagation

I'm confused by a couple of things in equation 9.

1. Shouldn't the numerator terms for the y**2 term be (k+1)? Not the 'k' index as shown?
2. Why is there a 2 in the third term - term for dU/dz? Isn't finite difference for dU/dz = (U(k+1) - U(k))/dz. The equation shows it as (U(k+1) - 2*U(k))/dz
3. I struggle to understand what A is in these equations. At first I thought it was index of refraction (e.g. ~= 1.0) However, it must be something different. Dimensionally, it must have units of 1/length**2.

Any chance there is a reference you could point towards for eqn. 8. Wikipedia is a good reference for Helmholtz equation in general, but I can't follow where the A*U term comes from.

[ldoyle]

Dear mikessut,

I'm afraid I can't help out much, the original code is basically as old as I am. However, during the porting to pure Python last year I also tried to understand what is going on, I hope my findings might point you in the right direction. (I will try to find my notes this week, remind me if I forget).

1. I think this is correct. The problem is one cannot solve for the next layer in both directions simultaneously, so we need to solve along X first (using the kth value in the other direction) and then repeating along Y (not shown as formula, here indices would be swapped, cf. explanation below formula). Otherwise the j-variation could not be collected away in factor f_i in eq. (10), making the differential equation harder/impossible to solve with the double sweep elimination method used.
2. This might indeed by a typo. If you compare with f_i in eq. (13), I believe there sould be no extra 2 in eq. (9).
3. I also struggled finding the exact form of A. I reconstructed it back from the actual code, I'll have to check my notes. I found a hint of what it should look like in an optics script in Google, keyword were something like "slowly varying envelope approximation" and "inhomogeneous media":
   https://www.iap.uni-jena.de/iapmedia/de/Lecture/Fundamentals+of+Modern+Optics1427752800/FoMO14_Script_2015_02_14s.pdf
   page 77
   However this is not what is calculated in the current LightPipes code, Fred could also not remember what the implemented approximations were exactly.

Regarding point 1., one may be able to do it better nowadays with more computing power. Especially since the double sweep implemented in pure Python is slower than the original Cpp version, I tried to find an alternative approach.
A true 2D-5-point method could work better, I found an example on scipy-cookbooks:
https://scipy-cookbook.readthedocs.io/items/discrete_bvp.html
However my first attempt failed miserably, which is why until now the code does exactly what has been used successfully in the past >20 years.

Do remind me if I forget to look for my notes.
Hope this helps, Lenny