Possible issue

[dibyendumajumdar]

Hi I am getting a bit strange result with this one:

%matplotlib inline
isdark = False
from rayoptics.environment import *
from rayoptics.elem.elements import Element
from rayoptics.raytr.trace import apply_paraxial_vignetting

# JP2019-152887 Example 3 (Nikon AF-S Nikkor 400mm f/2.8 E FL ED VR)
# Obtained via https://www.photonstophotos.net/GeneralTopics/Lenses/OpticalBench/OpticalBenchHub.htm

opm = OpticalModel()
sm  = opm.seq_model
osp = opm.optical_spec
pm = opm.parax_model
osp.pupil = PupilSpec(osp, key=['image', 'f/#'], value=2.89)
osp.field_of_view = FieldSpec(osp, key=['object', 'angle'], flds=[0., 3.14])
osp.spectral_region = WvlSpec([(486.1327, 0.5), (587.5618, 1.0), (656.2725, 0.5)], ref_wl=1)
opm.system_spec.title = 'JP2019-152887 Example 3 (Nikon AF-S Nikkor 400mm f/2.8 E FL ED VR)'
opm.system_spec.dimensions = 'MM'
opm.radius_mode = True
sm.gaps[0].thi=1e10
sm.add_surface([1200.37,5,1.5168,63.88])
sm.ifcs[sm.cur_surface].max_aperture = 69.61
sm.add_surface([1199.79,1])
sm.ifcs[sm.cur_surface].max_aperture = 69.61
sm.add_surface([207.079,17.5,1.43384,95.26])
sm.ifcs[sm.cur_surface].max_aperture = 68.61
sm.add_surface([-1127.53,44.9])
sm.ifcs[sm.cur_surface].max_aperture = 68.61
sm.add_surface([175.97,18,1.43384,95.26])
sm.ifcs[sm.cur_surface].max_aperture = 57.61
sm.add_surface([-397.271,3.07])
sm.ifcs[sm.cur_surface].max_aperture = 57.61
sm.add_surface([-360.24,6,1.61266,44.46])
sm.ifcs[sm.cur_surface].max_aperture = 55.36
sm.add_surface([353.184,90])
sm.ifcs[sm.cur_surface].max_aperture = 55.36
sm.add_surface([66.4844,4,1.795,45.32])
sm.ifcs[sm.cur_surface].max_aperture = 34.61
sm.add_surface([45.9182,15,1.49782,82.54])
sm.ifcs[sm.cur_surface].max_aperture = 31.61
sm.add_surface([1114.11,18.503])
sm.ifcs[sm.cur_surface].max_aperture = 31.61
sm.add_surface([2992.55,2.5,1.755,52.34])
sm.ifcs[sm.cur_surface].max_aperture = 24.11
sm.add_surface([118.04,3.35])
sm.ifcs[sm.cur_surface].max_aperture = 23.36
sm.add_surface([-241.694,3.5,1.84668,23.83])
sm.ifcs[sm.cur_surface].max_aperture = 23.36
sm.add_surface([-86.4136,2.4,1.53996,59.52])
sm.ifcs[sm.cur_surface].max_aperture = 23.36
sm.add_surface([64.2643,38.179])
sm.ifcs[sm.cur_surface].max_aperture = 23.36
sm.add_surface([0,1.5])
sm.set_stop()
sm.ifcs[sm.cur_surface].max_aperture = 18.979
sm.add_surface([90.0336,7.6,1.48749,70.43])
sm.ifcs[sm.cur_surface].max_aperture = 19.36
sm.add_surface([-63.8039,1.2])
sm.ifcs[sm.cur_surface].max_aperture = 19.36
sm.add_surface([-65.9768,1.9,1.84668,23.83])
sm.ifcs[sm.cur_surface].max_aperture = 18.36
sm.add_surface([-114.876,5])
sm.ifcs[sm.cur_surface].max_aperture = 18.36
sm.add_surface([300.359,3.5,1.84668,23.83])
sm.ifcs[sm.cur_surface].max_aperture = 18.36
sm.add_surface([-128.056,1.9,1.59319,67.94])
sm.ifcs[sm.cur_surface].max_aperture = 18.36
sm.add_surface([53.9004,3.1])
sm.ifcs[sm.cur_surface].max_aperture = 18.36
sm.add_surface([-347.542,1.9,1.755,52.33])
sm.ifcs[sm.cur_surface].max_aperture = 17.36
sm.add_surface([94.5337,4.19])
sm.ifcs[sm.cur_surface].max_aperture = 17.36
sm.add_surface([118.353,3.5,1.7725,49.68])
sm.ifcs[sm.cur_surface].max_aperture = 16.86
sm.add_surface([-384.382,0.1])
sm.ifcs[sm.cur_surface].max_aperture = 16.86
sm.add_surface([67.4622,4.5,1.64,60.14])
sm.ifcs[sm.cur_surface].max_aperture = 17.11
sm.add_surface([-340.421,1.9,1.84668,23.83])
sm.ifcs[sm.cur_surface].max_aperture = 17.11
sm.add_surface([246.642,6.5])
sm.ifcs[sm.cur_surface].max_aperture = 17.11
sm.add_surface([0,1.5,1.5168,63.88])
sm.ifcs[sm.cur_surface].max_aperture = 17.86
sm.add_surface([0,74.22])
sm.ifcs[sm.cur_surface].max_aperture = 17.86
sm.list_surfaces()
sm.list_gaps()
sm.do_apertures = False
opm.update_model()
apply_paraxial_vignetting(opm)
layout_plt = plt.figure(FigureClass=InteractiveLayout, opt_model=opm, do_draw_rays=True, do_paraxial_layout=False,
                        is_dark=isdark).plot()
sm.list_model()
# List the optical specifications
pm.first_order_data()
# List the paraxial model
pm.list_lens()
# Plot the transverse ray aberrations
abr_plt = plt.figure(FigureClass=RayFanFigure, opt_model=opm,
          data_type='Ray', scale_type=Fit.All_Same, is_dark=isdark).plot()
# Plot the wavefront aberration
wav_plt = plt.figure(FigureClass=RayFanFigure, opt_model=opm,
          data_type='OPD', scale_type=Fit.All_Same, is_dark=isdark).plot()
# Plot spot diagrams
spot_plt = plt.figure(FigureClass=SpotDiagramFigure, opt_model=opm,
                      scale_type=Fit.User_Scale, user_scale_value=0.1, is_dark=isdark).plot()

[mjhoptics]

The only odd thing I see is that the aberration in the blue is far out of focus. I believe this is because several of the glasses used in this design are ED glasses, i.e. anomalous dispersion glasses. For instance, 433.952 probably corresponds to Hoya FDC100 or Ohara S-FPL53. These are anomalous dispersion glasses and are not modeled well by just specifying nd and Vd. You need to know the partial dispersion as well and this is not normally specified in a patent design. The index/v-number specification is called a fictitious glass by CODE V, and they assume the fictitious glass corresponds to a "normal" partial dispersion. I do something similar but for these glass types, it's a poor approximation. I've attached a .roa file that has the 3 lenses with anomalous dispersion replaced by nearby Hoya glasses. I did the replacement modifying the curvatures to maintain the first order properties, since the glasses aren't an exact match. The aberrations are much better after doing this. (Strip the .txt off of the filename, that was the only way I could attach it). issue25.roa.txt

[dibyendumajumdar]

I was wondering if I would get better results using the glass types - there are fluorite elements and ED glasses in that lens.

If I wanted to use fluorite or Hikari glass catalogue (this being Nikon) - how can I do that?

p.s. the fluorite elements are described here: https://www.nikon.com/products/components/assets/pdf/caf2-e.pdf

[dibyendumajumdar]

Maybe there are similar issues in the other lenses too - for recent lenses I can probably find the correct matching glass - but for older lenses this might be a problem.

[mjhoptics]

I would use the InterpolatedGlass class in rayoptics.seq.medium. This takes pairs of wavelengths and index and uses scipy to interpolate them.

[dibyendumajumdar]

Is the issue restricted to some special glass types or is this a general issue?

I am looking at specifying more specific glass types where I can find an exact match from patent to glass maker, but I find that I cannot find a match for every element, and even if the glass has same abbe number and index as the spec, the output changes. That seems to indicate that the approximation used from abbe number/refractive index is a bit far off?

[mjhoptics]

yes, the problem will be more severe for ED glasses, but for high performance optics, replacing a fictitious glass with a real one usually required a re-optimization to regain the previous performance. For high performance systems we'd do melt adjustments, i.e. slightly re-optimizing a system (perhaps only changing the airspaces) to compensate for index changes due to the actual melt of glass used by the fabricator. Optical specifications can be quite sensitive to to index or dispersion variation; the tolerances on index can be as tight as 1e-4.

[dibyendumajumdar]

I opened a PR for a lens I previously submitted. I added another spec where I found the Hoya glass types used in the patent.

https://github.com/mjhoptics/ray-optics-notebooks/pull/5

It seems that the blue aberration is pretty off when I use just refractive index and abbe number. Using the Hoya glass types made a big difference.

Unfortunately I cannot always map to a glass type.

I will check if I can improve the Nikkor 400mm example by using InterpolatedGlass as you suggest.

[mjhoptics]

You could try using the glassmap application that is part of my opticalglass package. I just updated it (on pypi, conda-forge lags by several hours) to version 0.7.3 with a change that displays the partial dispersions of the 6 supported catalogs from F to d, i.e. the partials in the blue. You can see in the attached screenshot the map of the partials vs. v-number; the lower left hand corner is where all of the ED glasses live. The so-called "normal line" would go through the majority of the glasses from lower center to upper right. This is how glasses specified only by index and v-number will be modeled, as conforming to the normal glass line.

[dibyendumajumdar]

Thank you - I will look at that. Is there a way to search for a glass by refractive index / abbe number ? (I mean in ray-optics scripting)

I assume that this will still not cover Calcium Fluoride or Silica glass?

[dibyendumajumdar]

How do I use the glass type. I got as far as this:

from rayoptics.seq.medium import InterpolatedGlass

flglass = InterpolatedGlass('CaF2', [486.1327, 1.43700, 587.5618, 1.43384, 656.2725, 1.43245])

But I am not sure how I can specify this in add_surface().

[mjhoptics]

Yes, you can search the catalogs by index and v-number. You use the opticalglass package, which is installed by rayoptics as a dependency. Here's how to get the index and v-number (plus some other stuff) for, e.g. the Schott catalog. from opticalglass.glassfactory import get_glass_catalog from opticalglass.glass import get_glass_map_arrays schott = get_glass_catalog('Schott') nd, vd, PFd, bc1, bc2, names = get_glass_map_arrays(schott, 'd', 'F', 'C') The return values are the central index, v-number, partial dispersion from F to d, coefficients to a quadratic model of the index variation, and the name. The currently available catalogs are: "CDGM", "Hikari", "Hoya", "Ohara", "Schott", "Sumita" I'd like to have a catalog mechanism, kind of like the notebook archive idea, where additional material definitions could be made available and maintained.

[mjhoptics]

Sorry, I thought I had pressed "comment" on the aboce comment last night but apparently not...

add_surface doesn't currently accept a material-ish instance. I'd use the index/v-number entry as before and then add the following: sm.gaps[sm.cur_surface].medium = flglass

[dibyendumajumdar]

Am I making a mistake here:

from rayoptics.seq.medium import InterpolatedGlass

flglass = InterpolatedGlass('CaF2', [486.1327, 1.43700, 587.5618, 1.43384, 656.2725, 1.43245])

I also tried:

from rayoptics.seq.medium import InterpolatedGlass

flglass = InterpolatedGlass('CaF2', [(486.1327, 1.43700), (587.5618, 1.43384), (656.2725, 1.43245)])

I get an error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-251073314e1b> in <module>
      1 from rayoptics.seq.medium import InterpolatedGlass
      2 
----> 3 flglass = InterpolatedGlass('CaF2', [(486.1327, 1.43700), (587.5618, 1.43384), (656.2725, 1.43245)])
      4 
      5 opm = OpticalModel()

C:\Software\Python37\lib\site-packages\rayoptics-0.6.5.post0.dev1+g3a0af96-py3.7.egg\rayoptics\seq\medium.py in __init__(self, label, pairs, rndx, wvls, cat)
    140             self.wvls = wvls
    141             self.rndx = rndx
--> 142         self.update()
    143 
    144     def __repr__(self):

C:\Software\Python37\lib\site-packages\rayoptics-0.6.5.post0.dev1+g3a0af96-py3.7.egg\rayoptics\seq\medium.py in update(self)
    159     def update(self):
    160         self.rindex_interp = interp1d(self.wvls, self.rndx, kind='cubic',
--> 161                                       assume_sorted=False)
    162 
    163     def glass_code(self):

C:\Software\Python37\lib\site-packages\scipy\interpolate\interpolate.py in __init__(***failed resolving arguments***)
    531 
    532             self._spline = make_interp_spline(xx, yy, k=order,
--> 533                                               check_finite=False)
    534             if rewrite_nan:
    535                 self._call = self.__class__._call_nan_spline

C:\Software\Python37\lib\site-packages\scipy\interpolate\_bsplines.py in make_interp_spline(x, y, k, t, bc_type, axis, check_finite)
    825     if nt - n != nleft + nright:
    826         raise ValueError("The number of derivatives at boundaries does not "
--> 827                          "match: expected %s, got %s+%s" % (nt-n, nleft, nright))
    828 
    829     # set up the LHS: the collocation matrix + derivatives at boundaries

ValueError: The number of derivatives at boundaries does not match: expected 1, got 0+0

[mjhoptics]

I think you need more data points. I use cubic spline interpolation and my guess is that 3 data points isn't sufficient for this. I'd use all the data you can find from the i-line (365nm) to the near IR (e.g. r-line at 706.5nm).

[dibyendumajumdar]

Adding another value helped. I tried with just the CaF2 elements - results are better than before but still the blue dots are off. It seems I need to use glass specs for all elements, but sadly Nikon lenses always tend to have 1-2 elements that use glasses not listed anywhere. I don't know if they just have internal glass types or it is a deliberate obfuscation.

[mjhoptics]

I'll admit most of the time I examined lens patents, I was using them as a starting point for a slightly different problem and always used an optimizer on them. Patents rarely patent the exact prescription they are using, the claims usually put ranges on different parameters to codify the invention. The patents are written by attorneys who aren't technical experts and they can have an obfuscating effect on the disclosure, even though it might not be intentional.

[dibyendumajumdar]

It doesn't help also that the patent has errors.

I tried my best to match with appropriate Hikari glass types - the vd was incorrect a few times.

https://github.com/mjhoptics/ray-optics-notebooks/pull/6

[dibyendumajumdar]

When using interpolated glass type, I guess that the indices provided for the wavelengths are used - i.e. no need to interpolate? Is the vd parameter at all needed if we can provide the refractive index per wavelength?

[dibyendumajumdar]

btw when you mention optimization - is that using ray-optics or some other software?

[mjhoptics]

For the InterpolatedGlass type, as well as the other catalog glass types, the rindex() function always uses the interpolating polynomial, I don't check to see if the requested wavelength matches the defining set of wavelengths. The v-number parameter isn't needed in this case. I used CODE V for many years as a designer, that's the optimization I'm referring to.

[dibyendumajumdar]

Okay thank you