quartiq / rayopt

Python optics and lens design, raytracing

RayOpt

Introduction

Optics design (lenses, cavities, gaussian optics, lasers). Geometric, paraxial, and gaussian ray tracing.

Installation

Install like any usual Python package using pip, easy_install, or plain setup.py. Anaconda packages from three operating systems and three current Python versions are available through Anaconda. Install with:

conda install -c https://conda.anaconda.org/jordens/channel/ci rayopt

The distribution already contains all materials from http://refractiveindex.info/.

More materials

More materials and lenses catalogs can be obtained from the freely available versions of Oslo and Zemax, copied to catalog/ and then parsed using rayopt/library.py.

Zemax

More materials and lenses catalogs can be obtained from the freely available versions of Oslo and Zemax, copied to catalog/ and then parsed using rayopt/library.py (see there for details on where the files are expected)

Get Zemax optics studio. You can either install the software or unpack it with innoextract. Depending on your chosen method the paths have to be adapted:

$ python -m rayopt.library \
Users/$USER/My\ Documents/Zemax/Glasscat \
Users/$USER/My\ Documents/Zemax/Stockcat

OSLO

For OSLO, download and install OSLO.:

get and install http://www.lambdares.com/images/OSLO/OSLO662_EDU_Installer.exe
$ python -m rayopt.library \
Users/Public/Documents/OSLO66\ EDU/bin/lmo \
Users/Public/Documents/OSLO66\ EDU/bin/glc

Examples

Usage examples are at in their own repository as IPython Notebooks, specifically also the Tutorial.

Notes

Distance

The choice of prescription specification is a little different from most other lens design and ray tracing programs. RayOpt associates with an element (surface):

• distance (or directional offset, measured in the global, unrotated coordinate system) of the element's apex relative to the previous element's apex.
• orientation (x-y-z Euler angles in the rotating frame) with respect to the directional offset
• element properties (type of element, curvature, aspheric and conic coefficients, focal length of an ideal element)
• optionally, the material after the element (behind the surface)

Ray data are given at (ray intercepts) or just after (direction cosines, paraxial slopes) the respective element unless stated otherwise (e.g. incidence angles).

The choice of associating the "distance to" and not the "thickness after" with a surface has several advantages: shifts, offsets, tolerances can be implemented in a more straight forward manner, ray propagation becomes more natural and efficient (transfer, intercept, refraction), ray data at the surfaces' apex planes does not need to be tracked. The "thickness after" does not have much meaning in ray trace data as it can only be used later when tracing toward the next element and its direction is typically ill defined. Compared to most other programs the distance data is the thickness data shifted by one element towards the object.

Object and Image

Object and image are located at the first (index 0) and last (index -1) surface respectively. This naturally allows tracking their positions, material and shape data and supports curved objects and images naturally. Further data like pupils data are maintained in the two Conjugate s of the System.

Therefore, a minimal system of a single lens consists of fours surfaces: object, the two lens surfaces (one of which can be the aperture stop) and the image surface. The offset data of the first (object) surface does play a role in ray tracing but it can be useful as it locates the global coordinate system's origin. The material of the last (image) surface is used as it can cause incident rays to be evanescent at the image surface. This can also be compared to other programs where the thickness of the image surface is never relevant or the material in object space and the position of the lens has to be tracked somewhere else depending on the implementation.

Literature

• Warren J. Smith: Modern Optical Engineering, McGraw-Hill 2000: concise and methods derivation from paraxial all the way to arbitrary ray tracing, with terminology explained and examples given
• Michael Bass (ed): Mandbook of Optics I and II, OSA/McGraw-Hill 1995: physical foundations, broad on optics, comprehensive on theory, some info on numerics, some example designs
• Daniel Malacara: Handbook of Optical Design, Marcel Dekker Inc. 1994: Introduction, Aberations, Examples, more info on terminology, especially in ray tracing programs and codes
• Daniel Malacara: Geometrical and Instrumental Optics, Academic Press Inc. 1988: less info about algorithms and numerical methds, more examples and use cases, speciality lens designs
• Robert R. Shannon: The Art and Science of Optical Design, Cambridge University Press 1997: many examples with Oslo and Zemax, not very thorough on numerical methods and foundations, good for material comparison with own codes.
• Michael J. Kidger: Intermediate Optical Design, SPIE Press 2004: info on optimization techniques and algorithms, higher order aberrations, lots of example designs
• Milton Laikin: Lens Design, CRC Press 2007: little bit of basic theory, lots of basic and paradigmatic example designs
• Oslo Optics manual and reference
• Zemax manual

rayopt's Issues

Update object angle/CRA according to image size

Is this already possible or even desired? For me it would be helpful for designing a lens with a certain FOV.

Question: how to do close focus scenario

Hi @jordens ,
I am trying to add a close focus scenario for the Z Noct Nikkor 58mm f0.98 but I am getting an error.
I am guessing that I need to tell rayopt that the object is not an infinity but I don't know how to do this.

https://nbviewer.jupyter.org/github/cameragossip/cameragossip.github.io/blob/master/lensdesigns/Nikkor-Z-Noct-58mm-f0.95.ipynb

Thanks and Regards
Dibyendu

Formula for aspherical surfaces

Hi,
I believe the aspheric surface is computed in Spheroid.surface_sag() method?
I would like to know the formula used. In particular:

r2 = np.einsum("...i,...i", xy, xy)
if self.curvature:
     c, k = self.curvature, self.conic
      e -= c*r2/(1 + np.sqrt(1 - (1 + k)*c**2*r2))

In this what does r2 represent? Is it the height squared?
Is c a reciprocal of the input radius?

In the case of the 58mm Noct that I have looking at the patent specifies this formula:

x = (h**2 / r) / [1+ {1− (1 + κ) · (h / r)**2 } 1/2 ]
+ A4h**4 + A6h**6 + A8h**8 + A10h**10 + A12h**12 + A14h**14
Here, h is the height in the direction perpendicular to the optical axis, 
x is the distance along the optical axis direction from the tangent 
plane of the apex of the aspheric surface at the height h to the 
aspheric surface, 
and κ is the conic constant. 

Comparing the two, I think that if h**2 = r2, and c = 1/r then the two formulas match?

Question: Which value should match Fno on a patent

We get output such as:

front, back f number: [1.45 1.2221]
front, back working f number: [ inf 1.5338]

Which of these values should match the Fno conventionally specified in patent literature?

Thanks and Regards

Third order aberrations not being properly computed

The conic constant is ignored during the calculation of the third order aberrations

bayesian optimization

https://github.com/fmfn/BayesianOptimization
http://nbviewer.jupyter.org/github/fmfn/BayesianOptimization/blob/master/examples/exploitation%20vs%20exploration.ipynb

https://sheffieldml.github.io/GPyOpt/

Diffractive, periodic or parallel optical elements?

Hi ,

thanks for this great package. I have one questions: Do you see any way to trace rays through a periodic structure like a DMD, SLM or any kind of grating ? Additionally, I had problems setting even two macro mirror side by side with help of offset and angle, which reflect part of the light in different directions.

Thanks
Max

No materials data is installed

I used the conda method to install. However, I cannot run any of the examples as the material data does not seem to be installed. Or is it a path issue?

Question: Impact of cover glass in front of a digital sensor

I mentioned before that with some patents I am not getting a good match on the f number. I think that pattern I am seeing is that when a design has a cover glass in front of the image plane (in digital cameras newer designs take into account the impact of the cover glass) then I see the discrepancy. So presumably rayopt is considering the cover glass as just a regular optical element - and I guess that this impacting the result somehow? Any thoughts?

I think that the performance figures are adversely affected when the F no doesn't get close.

I have more examples you can see here : https://nbviewer.jupyter.org/github/cameragossip/cameragossip.github.io/tree/master/lensdesigns/

double_gauss example

Getting 'material None/None/sk4 not found'

GaussianTrace fails on "folded" systems

Trying to do a GaussianTrace on a system with folding mirrors. It will work if I don't fold the system, but the GaussianTrace plotting fails when I have the fold mirrors in the prescription.

It looks like elements' rot_axis isn't being set in some cases. If the element is straight but not normal, rot_axis won't be set, but is used in the from_axis method, which is used at least in GaussianTrace.

The below returns

AttributeError: 'Spheroid' object has no attribute 'rot_axis'

MWE:

import rayopt as ro
import matplotlib.pyplot as plt
import numpy as np

plt.close('all')

src = ro.Spheroid(material = ro.vacuum,
                  radius = 10.0e-3,
                  distance = 0.0)


m1 = ro.Spheroid(material = ro.mirror,
                 radius = 20.0e-3,
                 offset = (0.0, 0.0, 0.1),
                 angles = (-np.pi/4.0, 0.0, 0.0))

m2 = ro.Spheroid(material = ro.mirror,
                 radius = 20.0e-3,
                 offset = (0.0, 0.1, 0.0),
                 angles = (np.pi/4.0, 0.0, 0.0))

det = ro.Spheroid(material = ro.vacuum,
                  radius = 10.0e-3,
                  offset = (0.0, 0.0, 0.1),
                  angles = (0.0, 0.0, 0.0))

optics = ro.System(wavelengths = [355.0e-9],
                   scale = 1e0)

optics.extend([src, m1, m2, det])


optics.object = ro.FiniteConjugate(radius = 3.0e-3,
                                   pupil = dict(type = 'na',
                                                na = 0.1))

optics.update()

fig, ax = plt.subplots()
optics.plot(ax)

#for idx, el in enumerate(optics):
#    print('Element {}:'.format(idx))
#    print('  rotated:  {}'.format(el.rotated))
#    print('  straight: {}'.format(el.straight))
#    print('  normal:   {}'.format(el.normal))
#    print('  _distance:   {}'.format(el._distance))
#    print('  _direction:   {}'.format(el._direction))
#    print('  _offset:   {}'.format(el._offset))
#    print('  _angles:   {}'.format(el._angles))
#    print('  rot_axis:   {}'.format(el.rot_axis))
    
# Straight is if the direction is (0,0,1)
# Normal is if angles are all zero
# Rotated is if it is either not straight or not normal

gt = ro.GaussianTrace(optics)

fig, ax = plt.subplots()
gt.plot(ax)
optics.plot(ax)

automatically propose stop

In paraxial_trace you write:

# TODO
    # * introduce aperture at argmax(abs(y_axial)/radius)
    #   or at argmin(abs(u_axial))

could this be done in an instance method e.g. paraxial_ray.find_stop or .find_aperture ?

Import LightTools libs

Hello,

Do you have an exemple to import Glass libraries of LightTools software?

Thanks

Optimizing a Doublet: unclear how some parameters work

Hi all! First, thanks for creating RayOpt, it looks like a wonderful library! I'm looking forward to playing around with it.

I had a few questions around how the optimizer works. Full disclosure, I'm a relative newbie at optics in general, so some of my confusion might be stemming from that :)

I grabbed the optimization code from the triplet example and tried to simplify it down to a long, air-spaced doublet. After playing around with parameters for a while, I got it to a point where it "works" (ok'ish spot diagram, etc)... but I'm not sure what all the various components actually do so it's hard to say how well it is working :) I'm happy to polish this up into an example optimizer demo for the notebook repo, complete with comments

My full script is here, and a few questions about how it works below:

variables.extend(PathVariable(s, (i, "curvature"), (-1/10, 1/10))
                 for i in (1, 2, 3, 4))

How do the "bounds" work here? The docs for scipy.minimize suggest this is the extent the variables can reach, but in practice this seems to be something closer to step size during optimization? E.g. the ROC of a surface is not limited to [-0.1,0.1] when using the above parameters

def get(s):
    s.update()
    s.paraxial.focal_length_solve(500)
    s.paraxial.propagate(start=-2)
    s.paraxial.refocus()
    s.paraxial.propagate(start=-1)
    s.paraxial.resize()
    s.paraxial.aberrations()

In the get() function, what is the purpose of s.paraxial.propagate(start=-2) and s.paraxial.propagate(start=-1)? Skimming the source, it seems these might be the surfaces to start propogating from... but I don't understand why negative numbers are being used? I suspect I'm entirely incorrect about what these are for.

More generally, is there a high-level explanation of why the sequences is focal-length -> prop -> refocous -> prop -> resize -> aberr ?

operands = [
    FuncOp(s, get),
    #FuncOp(s, lambda s: s[-1].offset[2:], min=60),
    #FuncOp(s, lambda s: s.edge_thickness()[1:], min=2),
    #FuncOp(s, lambda s: np.diff(s.track), min=2),
    FuncOp(s, lambda s: s.paraxial.transverse3[:, 5:].sum(0)/.5, min=-1, max=1),
    FuncOp(s, lambda s: s.paraxial.transverse3[:, 3:5].sum(0)/1, min=-1, max=1),
    FuncOp(s, lambda s: s.paraxial.transverse3[:, :3].sum(0)/.03, min=-1, max=1),
    #FuncOp(s, lambda s: s.paraxial.transverse3[:, (0, 1, 2, 5, 6)].sum(0), weight=1),
    GeomOp(s, rays=16, weight=1),
    #PolyOp(s, weight=1),
]

My understanding is that the operands define the function that we are trying to minimize. So the get() ensures we are hitting the correct focal length, the commented-out lines for offset and thickness can be used to control details about the surfaces, the s.paraxial.transverse3 conditions have something to do with minimizing aberrations and spot size, and not sure what the GeomOp does.

Is there an easy explanation what the s.paraxial.transverse3 conditions do? I tried to dig around in the code but wasn't able to figure it out. Similarly, I'm not sure what the commented out sum(0) condition does, or the PolyOp

Thanks for the help! As I said, happy to take this and help put together an "optimizer 101" example for other clueless newbies like myself :)

Saving all intersection points for e.g. the spot plots into the Analysis object?

Is it already possible to obtain all intersection points of a given ray with all surfaces after the tracing?
As far as I see the information is at least for the GeometricRay.propagate dropped. How would one best save them for later use. E.g. computing phase at a given point.

pupil aiming with tilted object

s1=ro.system_from_yaml('''!!python/unicode 'description': !!python/unicode 'AC508-200-B AC508-200-B NEAR IR
  ACHROMATS: Infinite Conjugate 200'
!!python/unicode 'elements':
- !!python/unicode 'angles': [0.08726646259971647, 0.0, 0.0]
  !!python/unicode 'material': basic/air
  !!python/unicode 'radius': 0.01401147303517
- {!!python/unicode 'curvature': 0.0019409937888199, !!python/unicode 'distance': 193.14,
  !!python/unicode 'material': SCHOTT/N-SF6HT, !!python/unicode 'radius': 25.4}
- {!!python/unicode 'curvature': 0.0091575091575092, !!python/unicode 'distance': 5.0,
  !!python/unicode 'material': SCHOTT/N-LAK22, !!python/unicode 'radius': 25.4}
- {!!python/unicode 'curvature': -0.0074626865671642, !!python/unicode 'distance': 8.2,
  !!python/unicode 'material': basic/air, !!python/unicode 'radius': 25.4}
- {!!python/unicode 'distance': 100.0, !!python/unicode 'material': basic/air}
- {!!python/unicode 'material': basic/air, !!python/unicode 'radius': 18.0}
- {!!python/unicode 'curvature': 0.0074626865671642, !!python/unicode 'distance': 100.0,
  !!python/unicode 'material': SCHOTT/N-LAK22, !!python/unicode 'radius': 25.4}
- {!!python/unicode 'curvature': -0.0091575091575092, !!python/unicode 'distance': 8.2,
  !!python/unicode 'material': SCHOTT/N-SF6HT, !!python/unicode 'radius': 25.4}
- {!!python/unicode 'curvature': -0.0019409937888199, !!python/unicode 'distance': 5.0,
  !!python/unicode 'material': basic/air, !!python/unicode 'radius': 25.4}
- !!python/unicode 'angles': [-0.08726646259971647, 0.0, 0.0]
  !!python/unicode 'distance': 193.33044385190678
  !!python/unicode 'material': basic/air
  !!python/unicode 'radius': 10.0
!!python/unicode 'image':
  !!python/unicode 'pupil': {!!python/unicode 'distance': -404.3669122413583, !!python/unicode 'radius': 36.37432010274401,
    !!python/unicode 'refractive_index': 1.0002750477973053, !!python/unicode 'update_radius': true}
  !!python/unicode 'radius': 0.01401147303517
  !!python/unicode 'type': !!python/unicode 'finite'
!!python/unicode 'object':
  !!python/unicode 'pupil': {!!python/unicode 'distance': 404.1764683894514, !!python/unicode 'radius': 36.374320102743994,
    !!python/unicode 'refractive_index': 1.0002750477973053, !!python/unicode 'update_radius': true}
  !!python/unicode 'radius': 10.0
  !!python/unicode 'type': !!python/unicode 'finite'
!!python/unicode 'pickups': []
!!python/unicode 'scale': 0.001
!!python/unicode 'solves': []
!!python/unicode 'stop': 5
!!python/unicode 'validators': []
!!python/unicode 'wavelengths': [8.0e-07]''')

ro.Analysis(s1)

Question: Why is paraxial free space propagation divided by refractive index?

def paraxial_matrix(self, n0, l):
        m = np.eye(4)
        # 4x4 block matrix, M = [[A, B], [C, D]], A is 2x2 tan sag
        d = self.distance

then:

        m[0, 2] = m[1, 3] = d/n0
        return n0, m

Isn't the normal propagation matrix only dependent on the distance?

Thanks
Max

Matplotlib deprecation of Axis `set_smart_bounds`

In Matplotlib method set_smart_bounds is deprecated and throws the following error when performing analysis

AttributeError: 'XAxis' object has no attribute 'set_smart_bounds'

Seems that this capability was deprecated in matplotlib==3.2 but throws the error with older versions as well.

https://matplotlib.org/3.2.2/api/_as_gen/matplotlib.axis.XAxis.set_smart_bounds.html

use quaternions and not euler angles

obviously

Implementation of Off-Axis Parabolic Mirrors

Hi,
thanks for providing this great library.

I'm trying to test and optimize an optical system that uses off-axis parabolic mirrors under 90° to relay an image. I did manage to get an imaging system with 90° reflection to work and also to produce the correct parabolic shape with conic constant=-1, but the optical axis is always at the pole of the parabola.

Is there any way to implement this component in rayopt correctly? And if not, how hard would it be to adapt the API to allow for this shape to be used?

Kind regards,

Artur

Assertion error in brent solver

Hi,
I am not sure why I am getting an assertion failure with this example:

https://nbviewer.jupyter.org/github/cameragossip/cameragossip.github.io/blob/master/lensdesigns/Nikkor-400-f2.8-FL-ED-JP2019-152887.ipynb

Also is an aperture stop mandatory? I have an example without aperture stop that is also failing similarly.

Note that there is no error if I limit the field to 0 i.e. on-axis.

Thanks and Regards
Dibyendu

Matplotlib deprecation warning and numpy TypeError

Hi,
working through the jupyter tutorial I came across this warning.

ro.Analysis(s)
/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py:163: MatplotlibDeprecationWarning: 
The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.
  ax.xaxis.set_smart_bounds(True)
/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py:164: MatplotlibDeprecationWarning: 
The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.
  ax.yaxis.set_smart_bounds(True)
/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py:163: MatplotlibDeprecationWarning: 
The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.
  ax.xaxis.set_smart_bounds(True)
/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py:164: MatplotlibDeprecationWarning: 
The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.
  ax.yaxis.set_smart_bounds(True)
/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py:163: MatplotlibDeprecationWarning: 
The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.
  ax.xaxis.set_smart_bounds(True)
/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py:164: MatplotlibDeprecationWarning: 
The set_smart_bounds function was deprecated in Matplotlib 3.2 and will be removed two minor releases later.
  ax.yaxis.set_smart_bounds(True)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-25-00344cd64651> in <module>
----> 1 ro.Analysis(s)

/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py in __init__(self, system, **kwargs)
     72             setattr(self, k, v)
     73         if self.run:
---> 74             self.run()
     75         if self.print:
     76             for t in self.text:

/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py in run(self)
    133                                    sharex=True, sharey=True, squeeze=False)
    134             self.figures.append(fig)
--> 135             self.spots(ax[::-1], self.system.fields)
    136 
    137         if self.plot_opds:

/usr/local/lib/python3.9/dist-packages/rayopt/analysis.py in spots(self, ax, heights, wavelengths, nrays, colors)
    271                 r = paraxial.airy_radius[1]/paraxial.wavelength*wi
    272                 t = GeometricTrace(self.system)
--> 273                 t.rays_point((0, hi), wi, nrays=nrays,
    274                              distribution="hexapolar", clip=True)
    275                 # plot transverse image plane hit pattern (ray spot)

/usr/local/lib/python3.9/dist-packages/rayopt/geometric_trace.py in rays_point(self, yo, wavelength, nrays, distribution, filter, stop, clip)
    208                    distribution="meridional", filter=None, stop=None,
    209                    clip=False):
--> 210         ref, yp, weight = pupil_distribution(distribution, nrays)
    211         self.rays(yo, yp, wavelength, filter=filter, stop=stop,
    212                   clip=clip, weight=weight, ref=ref)

/usr/local/lib/python3.9/dist-packages/rayopt/utils.py in pupil_distribution(distribution, nrays)
    185         l = [np.zeros((2, 1))]
    186         for i in np.arange(1, n + 1.):
--> 187             a = np.linspace(0, 2*np.pi, 6*i, endpoint=False)
    188             l.append([np.sin(a)*i/n, np.cos(a)*i/n])
    189         xy = np.concatenate(l, axis=1).T

<__array_function__ internals> in linspace(*args, **kwargs)

/usr/lib/python3/dist-packages/numpy/core/function_base.py in linspace(start, stop, num, endpoint, retstep, dtype, axis)
    111 
    112     """
--> 113     num = operator.index(num)
    114     if num < 0:
    115         raise ValueError("Number of samples, %s, must be non-negative." % num)

TypeError: 'numpy.float64' object cannot be interpreted as an integer

Thanks for sharing!

microscope example

AttributeError: 'FullTrace' object has no attribute 'rays_paraxial_point'

Old version of FullTrace?

Mirror is infinite, when its finite.

It seem like rays of a parallel ray bundle reflect at the continuation of a mirror surface, even if the radius is smaller than the bundle diameter?

Question: is there a way to scale a design?

Hi,

Sometimes in patent literature the data is presented in a different scale than the actual values.
For instance the patent JP2019090947A describes the Z Nikkor 35mm f1.8 but presents the data scaled so that it appears to be a miniature lens. i would like to normalize so that different designs are comparable.

Thank you

Regards

Error running tutorial

I've just installed rayopt from the Conda repository:
conda install -c https://conda.anaconda.org/jordens/channel/ci rayopt

This is what conda tells me I have:
rayopt * 0.0+git py27_102 jordens/channel/ci

My python, numpy, and matplotlib are conda's latest
numpy 1.11.1
matplotlib 1.5.1
python 2.7.12

And finally OS X 10.11.5

I ran the tutorial jupyter notebook that comes with the package, but an error was raised in the third notebook frame (ro.Analysis(s))

Any ideas?

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-17-bb2a1fe7e9c5> in <module>()
----> 1 ro.Analysis(s)

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/analysis.pyc in __init__(self, system, **kwargs)
     72             setattr(self, k, v)
     73         if self.run:
---> 74             self.run()
     75         if self.print:
     76             for t in self.text:

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/analysis.pyc in run(self)
    113         for h in min(self.system.fields), max(self.system.fields):
    114             t = GeometricTrace(self.system)
--> 115             t.rays_clipping((0, h))
    116             t.plot(ax)
    117 

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/geometric_trace.pyc in rays_clipping(self, yo, wavelength, axis)
    213 
    214     def rays_clipping(self, yo, wavelength=None, axis=1):
--> 215         z, p = self.system.pupil(yo, l=wavelength, stop=-1)
    216         yp = np.zeros((3, 2))
    217         yp[1:, axis] = p[:, axis]/np.fabs(p).max()

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/system.pyc in pupil(self, yo, l, stop, **kwargs)
    592             c = self._pupil_cache[k] = PolarCacheND(self._aim_pupil,
    593                                                     l=l, stop=stop, **kwargs)
--> 594         q = c(*yo)
    595         return q[0], q[1:].reshape(2, 2)

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/cachend.pyc in __call__(self, *args)
     50             if np.any(np.isnan(guess)):
     51                 guess = self.guess
---> 52         value = self.solver(*args, guess=guess, **self.kwargs)
     53         self.cache[args] = value
     54         self._update()

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/system.pyc in _aim_pupil(self, xo, yo, guess, **kwargs)
    577             yp = [0, 0]
    578             yp[ax] = 2*sig - 1.
--> 579             a1 = self.aim_marginal(y, yp, z, a[sig, ax], **kwargs)
    580             a[sig, ax] = a1
    581             if sig == 1:  # and guess is None

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/system.pyc in aim_marginal(self, yo, yp, z, p, l, stop, **kwargs)
    553             else:
    554                 return d[-1]
--> 555         a = self.solve_brentq(dist, **kwargs)
    556         assert a
    557         return a*p

/Users/gpajer/miniconda2/lib/python2.7/site-packages/rayopt/system.pyc in solve_brentq(self, merit, a, b, tol, maxiter)
    496                 break
    497         if i == maxiter - 1:
--> 498             raise ValueError("no viable interval found", a, b, fb)
    499         fa = merit(a)
    500         if abs(fa) <= tol:

ValueError: (u'no viable interval found', 4.0, 4.0, nan)

Library.get_all: material None/None/sk16 not found

Hi I like the rayopt library and the integration with the scipy.optimize for lens curvature optimization. However, when I run the tutorial Triplet.ipynb in the rayopt-notebook repo, it fails to retrieve the material Schott/SK4. What is missing in the installation?

Reference: https://github.com/quartiq/rayopt-notebooks/blob/master/triplet.ipynb

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-2-ab0625f0217a> in <module>
     55 """
     56 
---> 57 s = system_from_text(text, columns.split(),
     58     description=description)
     59 s.object.angle = np.deg2rad(20)

~/Projects/raytracing/venv/lib/python3.8/site-packages/rayopt/formats.py in system_from_text(text, *args, **kwargs)
     83     n = max(len(l) for l in array)
     84     array = [l for l in array if len(l) == n]
---> 85     return system_from_array(array, *args, **kwargs)
     86 
     87 

~/Projects/raytracing/venv/lib/python3.8/site-packages/rayopt/formats.py in system_from_array(data, columns, shifts, material_map, **kwargs)
     74             mat = try_get(line, columns, "material")
     75             mat = material_map.get(mat, mat)
---> 76             m = Material.make(mat)
     77             el.material = m
     78     return s

~/Projects/raytracing/venv/lib/python3.8/site-packages/rayopt/material.py in make(cls, name)
    116         from .library import Library
    117         lib = Library.one()
--> 118         return lib.get("material", name, catalog, source)
    119 
    120     def __str__(self):

~/Projects/raytracing/venv/lib/python3.8/site-packages/rayopt/library.py in get(self, *args, **kwargs)
    118 
    119     def get(self, *args, **kwargs):
--> 120         for k in self.get_all(*args, **kwargs):
    121             return k
    122 

~/Projects/raytracing/venv/lib/python3.8/site-packages/rayopt/library.py in get_all(self, typ, name, catalog, source, **kwargs)
    132         res = res.order_by(Typ.name)
    133         if not res.count():
--> 134             raise KeyError("{} {}/{}/{} not found".format(
    135                 typ, source, catalog, name))
    136         for item in res:

KeyError: 'material None/None/sk16 not found'

in util.pupil_distribution , catch exception if unknown distribution is specified

right now the code throws

UnboundLocalError: local variable 'xy' referenced before assignment

because all the if statements are skipped.

help with simulation

Update
Issue has been resolved, there was a communication issue with RMS; RMS geometric spot size or RMS OPD. This has been resolved. The model and Rayopt now seem to give similar results, see final comment.

Dear @jordens,

I really like Rayopt. At the moment, I am trying to test the library by verifying the properties of a bundle with a diameter of 1.2 mm which is focused with an achromatic doublet and subsequently refracted by a tilted plane parallel plate. I have verified the longitudinal displacement, transversal displacement, the spot size and Rayleigh length with an analytical calculation based upon Basis Wavefront Aberration Theory for Optical Metrology, Wyant.
I am, however, not sure on how to calculate the RMS value with Rayopt. What should I do with the filter, the stop or clip!? What should I use as reference in my RMS!? The RMS values are not the same? This can not be explained by assuming the unit is wavelength RMS and millimeter RMS.

My Rayopt simulation is as follows;

%pylab inline
import warnings
import numpy as np
import matplotlib.pyplot as plt
import rayopt as ro

# Lens used 12.5mm Dia. x 90mm FL, VIS-NIR, Inked, Achromatic Lens from Edmund Optics
# LINK: http://www.edmundoptics.com/document/download/391099
filename='zmax_49332ink.zmx'   
with open(filename) as file:
    data=file.read()

# Parameters:
wavelength=405e-9        # wavelength [m]
D=1.2                             # diameter bundle [mm] see, s.scale=0.001 [m]
T=35                              # plate thickness [mm]
utilt=np.radians(24)       # tilt angle plate [radians]

# Radius plate
#  I can't remember why I use tangent, should not matter 
#   as long diameter is large enough
spol=T*np.tan(np.pi/8)

# Create the system
s=ro.zemax.zmx_to_system(data)
s.object.pupil.radius = D/2
# Ensures rays created with function ray_point are in the [-D/2,D/2] range
s.object.pupil.update_radius = False 
s.object.angle = np.radians(0)   # [radians]                     
s.wavelengths = [wavelength]
s.update()  
# changes needed to make the Zemax data compatible with Rayopt
del s[0]
# set physical size of the offset surface, i.e. the left line in the drawing
s[0].radius = 20  # [mm]
# sets the length between the first virtual offset surface and the lens
s[1].distance = 0  # [mm]
# add parallel plate to the system
s.insert(4,ro.elements.Spheroid(distance=10,material='SCHOTT/N-BK7',
                                                   diameter=spol*2,angles=[utilt,0,0])) 
s.insert(5,ro.elements.Spheroid(distance=T/np.cos(utilt),material='basic/air',
                                                   diameter=spol*2,angles=[utilt,0,0]))
#NOTE:  due to rotation the thickness increases to T/np.cos(utilt) 
#             if this is not done the transversal focus shift displacement
#             does not agree with the theoretical model
s.update()
#s.align(s.paraxial.n) # used by jordens for astigmatic focus shift, destroys rotation
q = ro.GaussianTrace(s)
#q.refocus()
# astigmatic focus shift , can also be obtained from print(q) and looking at table
#print("Astigmatic focus shift "+str(abs(q.waist_position.T[0][-1])-abs(q.waist_position.T[1][-1])))+" mm.")
print("The Gaussian waist radius is "+str(round(q.spot_radius[-1][0]*1000,2))+" micrometers.")
print("The Rayleigh range is "+str(q.rayleigh_range[-1][0])+ " mm.")
# Geometric trace
g = ro.GeometricTrace(s)
# In my system, I am only interested in one field
# with a field angle equal to zero radians
# Several distribution can be chosen; hexapolar, random, radau
# The radau scheme should be able to give a good result while using not so many rays
fieldangle=0
g.rays_point((0, fieldangle), wavelength=wavelength, nrays=20, 
                     distribution="radau", filter=False, clip=False)
# Geometric focus [used]
g.refocus()
# The geometric RMS spotsize is then calculated at the focal point
# i.e. RMS= <(W-<W>)2>1/2
# on default Rayopt specifies the focal point at the last surface 
# as it sets i=surface equal to -1.
# all rays are given the same "weight"
print("RMS geometric spotsize is "+str(g.rms()*1000)+" micrometers.")
# The focus point distance is measured with respect to the lens
print("The focus point distance from the lens is "+str(g.path[-1]-g.path[3])+" mm.")
print("The transversal displacement is "+str(g.y[-1,-1,1])+" mm.")
p, q, opd = g.opd(resample=False)
print("The lambda OPD RMS is "+str(np.sqrt((opd**2 * g.w).sum()/g.w.sum())))
#
p = ro.ParaxialTrace(s)
print("The Airy radius is "+ str(p.airy_radius[1]*1000)+" micrometers.")
# paraxial focus [not used]
#s.paraxial.refocus()
ro.Analysis(s,refocus_full=False, update=False)
# Gaussian trace
#   plot only works at ultilt is 0 degrees
if utilt==0:  
    fig, ax = plt.subplots()
    s.plot(ax)
    q.plot(ax, color="red", scale=1)
# Seidel aberrations
#z = ro.PolyTrace(s)
#str(z)
# Retrieve seidel
#print("\n".join(z.print_seidel()))

My analytical calculations are as follows:

import numpy as np
from scipy.integrate import dblquad

# Parameters used:
n=1.53                       #  should equal s.refractive_index(405e-9,4), where s is system defined above
utilt=np.radians(24)   # tilt angle [radians]
wavelength=0.405     # wavelength [micrometers]
T=35                          # plate thickness [mm]

# Lens used 12.5mm Dia. x 90mm FL, VIS-NIR, Inked, Achromatic Lens from Edmund Optics
D=1.2                # diameter bundle [mm]
efl=90               # effective focal length [mm]
bfl=88.47            # back focal length [mm]
ftol=0.02            # focal length tolerance, i.e. ftol*100 percent 
d=12.5               # lens diameter [mm]

# F-number
# source: https://en.wikipedia.org/wiki/F-number
fnumber=efl/D
# Spot size
# source: https://www.newport.com/n/gaussian-beam-optics
waist=2*wavelength/3.14*fnumber
print("The spot radius is "+str(round(waist,2))+" micrometers.")
# Rayleigh range see https://en.wikipedia.org/wiki/Rayleigh_length 
raileigh=np.pi*waist**2/(wavelength*1000)
print("The Rayleigh range is "+str(round(raileigh,3))+" mm.")
# Longitudinal focus shift
#   Wyant page 41 equation 68
slong=(n-1)/n*T
print("The focus point is at "+str(slong+bfl)+' mm.')
# Transversal focus shift
#   Wyant page 41 equation 70
disp=T*np.sin(utilt)*(1-np.sqrt((1-np.sin(utilt)**2)/(n**2-np.sin(utilt)**2)))
print("The transversal displacement is "+str(disp)+" mm.")
# 3rd order Seidel aberrations in [mm]
#  Wyant page 42 equation 72; spherical aberration
sabr=-T/pow(fnumber,4)*((pow(n,2)-1)/(128*pow(n,3)))
#   Wyant page 44 equation 75; coma
#   Note: cosine has been omitted, will be added back when function f is defined
coma=-T*utilt/pow(fnumber,3)*((pow(n,2)-1)/(16*pow(n,3)))
#   Wyant page 45 equation 77 ; astigmatism
#   Note: cosine has been omitted, will be added back when function f is defined
astig=-T*pow(utilt,2)/pow(fnumber,2)*((pow(n,2)-1)/(8*pow(n,3)))
# To calculate the combined RMS we add back the functional forms to the coefficients
# They have been obtained from Wyant, page 17, table 2
# The function f defines the wavefront aberration denoted by Wyant as w
def f(theta,rho):
	return sabr*(rho**4)+astig*(rho**2)*(np.cos(theta)**2)+coma*(rho**3)*np.cos(theta)
# We now try to evaluate equation 62, page 37, Wyant
# First we define two auxiliary functions
def ws(theta,rho):
   return f(theta,rho)**2*rho
def w(theta,rho):
   return f(theta,rho)*rho
# Hereafter we evaluate the integral, i.e. equation 62
var=1/np.pi*dblquad(ws,0,1,lambda rho: 0, lambda rho: 2*np.pi)[0]\
-1/(np.pi**2)*(dblquad(w,0,1,lambda rho: 0, lambda rho: 2*np.pi)[0]**2)
# the RMS value is obtained by the root
rms=np.sqrt(var)
# The lambda RMS is expressed in [mm] and will be converted to wavelength units
lambdarms=rms/(wavelength*1E-3)
print("The lambda OPD RMS is " "+str(round(lambdarms,6)))
# The Strehl radius is calculated with the first three terms of its Taylor series
# Wyant page 39 equation 67
rstrehl=1-(2*np.pi*lambdarms)**2+(2*np.pi*lambdarms)**4/2
print("The Strehl radius is "+str(round(rstrehl,2)))

elphel_v2

ValueError: ('no viable interval found', 0.0, 9.313225746154785e-10, nan)

yaml deprecation warning

Hi,
working through the jupyter tutorial I came across this warning.

s = ro.system_from_yaml(...)
/usr/local/lib/python3.9/dist-packages/rayopt/formats.py:89: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
  dat = yaml.load(text)

It was easy to suppress it, but there might be a better way.

import yaml
yaml.warnings({'YAMLLoadWarning': False})

Thanks for sharing!

Installing rii materials?

I have trouble running the example files because none of the materials appear to have installed when I ran setup.py. Any suggestions?

unnecessary condition in paraxial trace?

https://github.com/jordens/rayopt/blob/master/rayopt/paraxial_trace.py#L62-L63

When would this condition hold, especially n[0].shape==n ?
Because one line before n is set to n=self.length

Thanks
Max

paraxial matrix for arbitrary tilts

def paraxial_matrix(self, n0, l):
        # Reflection and Refraction of Gaussian Light Beams at
        # Tilted Ellipsoidal Surfaces
        # G. A. Massey and A. E. Siegman
        # Applied Optics IP, vol. 8, Issue 5, p.975
        #
        # [y', u'] = M * [y, u]
        n, md = super(Spheroid, self).paraxial_matrix(n0, l)
        c = self.curvature
        if self.aspherics is not None:
            c = c + 2*self.aspherics[0]
        theta = self.angles[0] if self.angles is not None else 0.
        costheta = np.cos(theta)
        m = np.eye(4)
        if self.material is not None:
            if self.material.mirror:
                m[2, 0] = 2*c*costheta
                m[3, 1] = 2*c/costheta
            else:
                mu = n/n0
                p = np.sqrt(mu**2 + costheta**2 - 1)
                m[1, 1] = p/(mu*costheta)
                m[2, 0] = n0*c*(costheta - p)
                m[3, 1] = mu*m[2, 0]/(costheta*p)
                m[3, 3] = 1/m[1, 1]
        m = np.dot(m, md)

        if self.angles is not None:
            # FIXME angles is incomplete:
            # rotate to meridional/sagittal then compute total incidence
            # angle, matrix, then rotate back
            phi = self.angles[2]
            cphi, sphi = np.cos(phi), np.sin(phi)
            r1 = np.array([[cphi, -sphi], [sphi, -cphi]])
            r = np.eye(4)
            r[:2, :2] = r[2:, 2:] = r1
            m = np.dot(r, np.dot(m, r.T))

Maybe you could give me hint. It is unclear to me why the angle[0] would be in anyway special. Why is angle[1] not taken into account, or is this indeed what is meant by the FIXME comment: that we rotate the system until only one angle is nonzero?
But: Why is theta in the above function equal angle[0] and not something like u[0] +angle[0]

Thx
Max

pip install error

trying to install rayopt via pip on Windows, with Python 3.6.

keep coming up with this error:
ValueError: 'rayopt/simplex_accel.pyx' doesn't match any files

Any tips?