bryancole / raypier_optics

A raytracing toolkit for optical design
Other
58 stars 8 forks source link

Square root of negative value in field calculation #14

Open quincyhuynh opened 1 year ago

quincyhuynh commented 1 year ago

Hi Bryan,

Your library has been super useful for evaluating optical systems!

However, I noticed the following when performing the field calculations in ./core/cfields.pyx, in line 94 of sum_gaussian_modes: inv_root_area = sqrt(sqrt(A.imag*C.imag -(B.imag*B.imag))*(2.0/M_PI))

In some situations A.imag*C.imag < (B.imag*B.imag) and the result of inv_root_area is Nan and the whole field calculation yields a matrix of NaN.

Would the right fix be to take fabs(sqrt(A.imag*C.imag -(B.imag*B.imag))?

Also I am curious which equations you're using to calculate the E_field.

Thanks!

bryancole commented 1 year ago

Hello Quincy,

I'm happy to hear that you're finding Raypier useful.

The A,B & C coefficients describe the gaussian mode via the Gamma 2x2 matrix described here(an "ABCD" matrix):

https://www.researchgate.net/publication/225711573_Simple_formula_for_a_Gaussian_beam_with_general_astigmatism_in_a_homogeneous_medium

There is no D coefficient because the matrix is symmetric (i.e. D==B). The real parts of these coefficients describe the wavefront curvature, the imaginary part describes the spatial width of the modes. A & C are the diagonal elements. a non-astigmatic mode would have B==0.

My suspicion is that if you have |A.imagC.imag < (B.imagB.imag),| then this mode will be in violation of the paraxial approximation. Thus, although you can avoid the NaNs using the change that you propose, the result may still be unphysical. It is important to remember than Gaussian modes are "nearly collimated beams". The depiction in text books where a Gaussian modes comes down to a nice narrow waist is rather misleading. It be a good approximate (and they are only an approximate) solution to Maxwell's Equations, the beam waist must be >> wavelength.

Mentally, I can't visualise what a mode with B.imag^2 > A.imag*C.imag looks like. This suggests a highly astigmatic mode. It may be worth visualising such a mode to get a feel if it is physically realistic or not.

In summary, yes, feel free to try with the 'fabs' function to avoid the Nans, but check the field profile carefully to see if it is still a solution to Maxwell's equations. How you do this check, I'm not sure off the top of my head. I'd need to sit down and think on it for a while. I think the 'lagrange invariant' is probably what you want to look at but the literature was not terribly helpful in explaining what the legrange invariant means (and maybe I'm not smart enough to figure it out). GaussletCollection objects have a legrange_invariant property to make it easy to explore this but I've not used/understood it that well.

I'm not going to have any time myself to look into this for the next month or two but feel free to submit a pull-request it you want some changes incorporated.

Best regards, Bryan

On 27/01/2023 04:07, Quincy Huynh wrote:

Hi Bryan,

Your library has been super useful for evaluating optical systems!

However, I noticed the following when performing the field calculations in |./core/cfields.pyx|, in line 94 of |sum_gaussian_modes|: |inv_root_area = sqrt(sqrt(A.imagC.imag -(B.imagB.imag))*(2.0/M_PI))|

In some situations |A.imagC.imag < (B.imagB.imag)| and the result of |inv_root_area| is |Nan| and the whole field calculation yields a matrix of |NaN|.

Would the right fix be to take |fabs(sqrt(A.imagC.imag -(B.imagB.imag))|?

Also I am curious which equations you're using to calculate the |E_field|.

Thanks!

— Reply to this email directly, view it on GitHub https://github.com/bryancole/raypier_optics/issues/14, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABYLNVSLMOZMSOM4IK6U2DWUNCYJANCNFSM6AAAAAAUIIAW4A. You are receiving this because you are subscribed to this thread.Message ID: @.***>

quincyhuynh commented 1 year ago

Hi Bryan,

Thanks for getting back to me! I will look into the gaussian beam formula more carefully and try to see what is causing this case and if it makes any physical sense.

Cheers, Quincy