Wrong root is picked in intersect function for Spherical surface profile.

[pwarciszewski]

Hi,

I found a bug in the following function: https://github.com/mjhoptics/ray-optics/blob/f055ded7e39ccf1fca54a1c1ec2f53f66a1705b0/src/rayoptics/elem/profiles.py#L287 It seems to be picking the wrong quadratic equation root and as a result, breaking the whole raytracing process. However, the commented-out function: https://github.com/mjhoptics/ray-optics/blob/f055ded7e39ccf1fca54a1c1ec2f53f66a1705b0/src/rayoptics/elem/profiles.py#L290 can handle the problem.

This issue arises when a simple relay lens, consisting of two positive lenses, is used and their separation causes the effective focal length to change from positive to negative. When this occurs, along with an aperture stop located in the middle of the two lenses, mathematical errors can occur within the intersect function. The following code provides an example:

from rayoptics.environment import *

isdark = False

opm = OpticalModel()
sm  = opm['seq_model']
osp = opm['optical_spec']
pm = opm['parax_model']
em = opm['ele_model']
pt = opm['part_tree']
ar = opm['analysis_results']

osp['pupil'] = PupilSpec(osp, key=['object', 'NA'], value=0.01)
osp['fov'] = FieldSpec(osp, key=['object', 'height'], value=5, flds=[0., 0.707, 1.], is_relative=True)
osp['wvls'] = WvlSpec([(780.2, 1.0)], ref_wl=0)

opm.radius_mode = True

sm.gaps[0].thi=100

sm.add_surface([100, 3, 'N-BK7', "Schott"]) 
sm.add_surface([-100, 98]) # when thi changed to 97, works
sm.add_surface([0,98]) # when thi changed to 97, works
sm.set_stop()
sm.add_surface([100, 3.1, 'N-BK7', "Schott"])
sm.add_surface([-100, 2])

opm.update_model()
pm.first_order_data()

layout_plt = plt.figure(FigureClass=InteractiveLayout, opt_model=opm, is_dark=isdark).plot()
plt.show()

While I can understand that is a peculiar situation a similar bug occurs also when only half of the system is considered:

from rayoptics.environment import *

isdark = False

opm = OpticalModel()
sm  = opm['seq_model']
osp = opm['optical_spec']
pm = opm['parax_model']
em = opm['ele_model']
pt = opm['part_tree']
ar = opm['analysis_results']

osp['pupil'] = PupilSpec(osp, key=['object', 'NA'], value=0.01)
osp['fov'] = FieldSpec(osp, key=['object', 'height'], value=5, flds=[0., 0.707, 1.], is_relative=True)
osp['wvls'] = WvlSpec([(780.2, 1.0)], ref_wl=0)

opm.radius_mode = True

sm.gaps[0].thi=100

sm.add_surface([100, 3, 'N-BK7', "Schott"]) 
sm.add_surface([-100, 98]) # when thi changed to 97, works
sm.add_surface([0,0])
sm.set_stop()

opm.update_model()
pm.first_order_data()

layout_plt = plt.figure(FigureClass=InteractiveLayout, opt_model=opm, is_dark=isdark).plot()
plt.show()

From what I understand, the problem is related to the inputs of the function. I noticed, that the d parameter is getting the opposite sign when the bug occurs, while z_dir stays the same. As the other super().intersect does not rely on z_dir, the problem does not appear. The issue may be related to the position of the aperture stop since the bug does not appear when it is placed on another surface. Because of that, I would guess it may originate already here:

https://github.com/mjhoptics/ray-optics/blob/f055ded7e39ccf1fca54a1c1ec2f53f66a1705b0/src/rayoptics/parax/firstorder.py#L187

when already for the first surface enp_dist is miscalculated, or not handled as intended. However, this is only a suspicion and not a definitive statement.

[mjhoptics]

Hi @pwarciszewski, Thank you for the detailed report. Your system is interesting because it is telecentric in object and image space. I've been working on this general part of the code on the develop branch and this is a useful testcase. I'll keep you updated on progress. Thanks again, Mike

[pwarciszewski]

Hi @mjhoptics First of all, thank you for the great Python module and its continuous development. I look forward to hearing more on the matter.

Pawel

[mjhoptics]

I've updated a fix to the master branch and merged it out to develop.

The problem had to do with ray aiming. The model provided has the interesting attribute that the entrance pupil is behind the object, i.e. it is farther from the optics than the object surface. Since the system is approximately telecentric in object space, this is not too surprising, but it is a corner case. The program was calculating the starting direction for a ray by constructing a ray going from the object point to the desired point in the pupil. The case where the pupil was virtual, as in your model, wasn't accounted for. The result was a ray that was traced away from the optics, which resulted in a failed ray trace (all nans for every surface). The remedy was to test the generated ray direction against the propagation direction (z_dir in object space) and reverse the ray direction to match z_dir, if needed.

As a side note, in the comments it was noted that an airspace of 97 (vs 98) didn't see a problem. In this case, the smaller spacing to the stop was enough to move the entrance pupil to the system side of the object surface so that it was no longer virtual.

The change was only needed in the ray iteration code; it had already been introduced in the trace_base() function.

Your observation that the ray trace was trying to find the wrong intersection point was spot on: the ray-trace was tracing the ray in the direction away from the lenses. Of course it's going to miss the surfaces. The root cause was upstream of the immediate failure, and had been accounted for in the main ray trace function, trace_base(), but not in iterate_ray()

[pwarciszewski]

Hi, it looks fine now, thank you for the quick update. I am happy that I could help.