help with simulation

[hstarmans]

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)))

[jordens]

rayopt already calculates primary (third order) aberration coefficients for you. See the example notebooks. http://nbviewer.jupyter.org/github/jordens/rayopt-notebooks/blob/master/focuser.ipynb

 # T       SA3      CMA3      AST3      PTZ3      DIS3     TACHC      TCHC
 0 S         0         0         0         0         0         0         0
 1 S   -0.1048  -0.00013-1.613e-07-2.757e-07-5.421e-10         0         0
 2 S   -0.1424 0.0002215-3.444e-07-6.616e-08 6.384e-10        -0         0
 3 S  -0.06914-7.929e-05-9.092e-08-5.932e-07-7.845e-10         0         0
 4 S   -0.1812 0.0003056-5.154e-07 1.815e-07  5.63e-10        -0         0
 5 S    0.1125 0.0001328 1.568e-07-9.931e-07-9.875e-10         0         0
 6 S     0.383 4.653e-05-6.765e-07 7.711e-07  7.85e-10         0         0
 7 S         0        -0         0         0         0        -0         0
     -0.002114 0.0004971-1.632e-06-9.755e-07-3.276e-10         0         0

And it can calculate aberration coefficients to arbitrary (!) order, see e.g. the end of http://nbviewer.jupyter.org/github/jordens/rayopt-notebooks/blob/master/tutorial.ipynb

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.

If you want to do geometric ray propagation and then calculate RMS spot size from that, have a look at the optimization examples in http://nbviewer.jupyter.org/github/jordens/rayopt-notebooks/blob/master/triplet.ipynb This is it:

t = GeometricTrace(system)
t.rays_point((0, field), nrays=nrays, distribution="radau", clip=False, filter=False)
print(t.rms())

[hstarmans]

Thank you for your quick reply: I have seen in the example that Rayopt calculates primary (third order) aberration coefficients. I realized, I put the s.object.angle at 0, you then get nan and no (third order) aberration coefficients are printed. This has been changed by using np.radians(0.02).

 # T       SA3      CMA3      AST3      PTZ3      DIS3     TACHC      TCHC
 0 S         0         0         0         0         0         0         0
 1 S  -0.01301-4.054e-05-1.264e-07-1.933e-07-9.962e-10         0         0
 2 S     0.051-5.748e-05 6.479e-08 5.354e-08-1.334e-10        -0         0
 3 S  -0.03484 9.705e-05-2.704e-07-1.112e-07 1.063e-09        -0         0
 4 S  0.006966-3.491e-05  1.75e-07        -0 -9.19e-10        -0         0
 5 S -0.004083 2.047e-05-1.026e-07         0  1.12e-09        -0         0
 6 S         0        -0         0         0         0        -0         0
      0.006036-1.543e-05-2.595e-07 -2.51e-07 1.348e-10         0         0

The figure looks as follows. I don't get what I have to put for my s.object.angle. I thought it was just the angle the chief ray makes, which should be zero !?

[jordens]

A geometric trace is only "more accurate" in the sense that it calculates all orders. But it is still limited by sampling error, i.e. limited number of rays. But the code I quoted above should be very close to what you want.

Your trace is giving you a bunch of NaN because your object size is zero (you have a zero angular radius of the infinite object). In principle rayopt should be able to work with that: the only non-zero aberration coefficient is the spherical aberration. Even third order astigmatism is technically zero because it is field-dependent: the plate would never be considered. But in the calculations rayopt assumes finite object size in a couple places and you get NaN.

[hstarmans]

The code already in Rayopt is already close to what i want :-). The only problem I see so far is that

g.rays_point((0, fieldangle), wavelength=wavelength, nrays=10, 
                     distribution="meridional", filter=False, clip=False)

Creates 11 rays in the range -6.25 mm up to 6.25 mm whereas I thought I had defined a pupil radius of D/2=0.6 i.e. so the 11 rays should be in the range -0.6 mm to 0.6 mm.

[jordens]

It will update the entrance pupil size according to the size of the stop in object space. I don't know what your stop is but I guess it's your first or second surface. You can insert an actual stop (i.e. a zero thickness, air-to-air interface) to stop your system down to the D=1.2 you want. But remember that the actual stop can't be in contact with the object plane. Otherwise aiming will fail.

[jordens]

... You can also system.object.pupil.update_radius = False; system.object.pupil.radius = .6.

[hstarmans]

That does help. The rays are now within the right range. Rayopt concludes the RMS is 0.007 mm at 24 degrees. The analytical calculations conclude it is 1.28e-05 mm that is 0.03 lambda RMS. The number which results from Rayopt seems too large by a factor of 547. An RMS of 0.007 mm is also quite large for a wavelength of 0.405 micron. Do I need to add a weighting for the rays!? I thought the default weighting of Rayopt which is one, should be fine.

[jordens]

I don't have your the data in front of me. And I can't see what could be wrong. If you g.refocus() after g.rays_point(...), does it change g.rms()?

[hstarmans]

g.refocus does influence the result it goes from 0.007 to 0.006 mm

System: 49332INK ACHROMAT
Scale: 1.0 mm
Wavelengths: 405 nm
Fields: 0
Object:
 Semi-Angle: 0 deg
 Pupil:
   Pupil Distance: 0
   Refractive Index: 1.00028
   Radius: 0.6
Image:
 Radius: 0.328
 Update Radius: True
 Pupil:
   Pupil Distance: -74.5488
   Refractive Index: 1.00028
   Update Radius: True
   Radius: 5.69628
Stop: 1
Elements:
 # T   Distance   Rad Curv   Diameter          Material       n      nd      Vd
 0 S          0        inf         40         basic/air   1.000   1.000   89.30
 1 S          0       55.8       12.5      SCHOTT/N-BK7   1.530   1.517   64.17
 2 S       1.98     -40.17       12.5      SCHOTT/N-SF5   1.711   1.673   32.25
 3 S        1.6     -116.3       12.5         basic/air   1.000   1.000   89.30
 4 S         10        inf     28.995      SCHOTT/N-BK7   1.530   1.517   64.17
 5 S     38.312        inf     28.995         basic/air   1.000   1.000   89.30
 6 S     54.698        inf    0.65636         basic/air   1.000   1.000   89.30

lagrange: -0
track length: 51.892
object, image height: [ 0.  0.]
front, back focal length (from PP): [ nan  nan]
entry, exit pupil height: [ nan  nan]
entry, exit pupil distance: [ nan  nan]
front, back focal distance: [ nan  nan]
front, back principal distance: [ nan  nan]
front, back nodal distance: [ nan  nan]
front, back numerical aperture: [ 0.     0.007]
front, back f number: [ nan  nan]
front, back working f number: [    inf  68.165]
front, back airy radius: [   inf  0.034]
transverse, angular magnification: [ -0.  nan]

 # T      path         n   axial y  axial nu   chief y  chief nu
 0 S         0         1       0.6         0        -0         0
 1 S         0      1.53       0.6 -0.005698         0         0
 2 S      1.98     1.711    0.5926 -0.003032         0         0
 3 S      3.58         1    0.5898 -0.006635         0         0
 4 S     13.58      1.53    0.5524 -0.006288         0         0
 5 S     51.89         1    0.3384 -0.007337         0         0
 6 S     98.03         1 2.362e-17 -0.007337         0         0

 # T       SA3      CMA3      AST3      PTZ3      DIS3     TACHC      TCHC
 0 S         0         0         0         0         0         0         0
 1 S       nan       nan       nan        -0       nan       nan       nan
 2 S       nan       nan       nan         0       nan       nan       nan
 3 S       nan       nan       nan        -0       nan       nan       nan
 4 S       nan       nan       nan        -0       nan       nan       nan
 5 S       nan       nan       nan         0       nan       nan       nan
 6 S       nan       nan       nan         0       nan       nan       nan
           nan       nan       nan         0       nan       nan       nan

Rayopt concludes: The focus point is at 102.27 mm (24 degrees), 101.03 mm (0 degrees) The transversal displacement is 5.42 mm (24 degrees) RMS is 0.0057mm (24 degrees) The spot radius is 19.42 micrometer (can only be calculated at 0 degrees) The Rayleigh range is 2.93 mm (can only be calculated at 0 degrees) Numerical model concludes: The focus point is at 100.6 mm. The transversal displacement is 5.42 mm (24 degrees) The combined lambda RMS value is 0.032 The spot radius is 19.35 micrometers. The Rayleigh range is 2.90 mm.

Concluding; - there is still a mismatch in RMS, Rayopt states a large RMS

• I have to tried to evaluate the RMS at a different point with the Rayopt model e.g. 100.6 instead of 102.27 mm (but this increases RMS, as it should with a good model)

[jordens]

Between which numbers do you think there is a mismatch? I have trouble figuring out what exactly you mean by those numbers and what you are comparing.

[hstarmans]

I am looking at the RMS value. Numerical model RMS equals 0.032 lambda which is 1.3e-05 mm Rayopt RMS is 0.0057 mm These numbers don't match, according to Rayopt there is a severe aberration. I am pretty sure, that the numerical model is correct. The RMS value of Rayopt seems too large. At 4 degrees tilt, Rayopt already has an RMS of 0.00014 mm; which is huge!

[jordens]

With your system I am seeing:

• 0.14 lambda peak-to-peak OPD
• 4 µm RMS geometric (!) spot size
• 37 µm Airy radius
• 19 µm Gaussian waist radius
• 2.5mm astigmatic focus shift (tangential waist to sagittal waist, longitudinally)

https://nbviewer.jupyter.org/github/jordens/sandbox/blob/master/rayopt-issue-16.ipynb

[jordens]

And this is with focusing for minimum geometric RMS. Your analytic formulae seem to do it for the paraxial focus. I am pretty sure that your 13 nm RMS spot size are incorrect. But I'd be happy to fix it if it is a bug in rayopt.

[hstarmans]

@jordens, thank you for your reply and model!! I will need some time to study the answer. I am looking at the simulation to understand if I can scan with a laser bundle if I tilt a transparent plate; this could be useful in a laser scanner.

[hstarmans]

I have looked at the system and improved my system with your comments, thanks! I am interested in the RMS optical path difference . Rayopt gives me the RMS geometric spot size. How do I calculate the RMS optical path difference with Rayopt!? The numbers don't match because they refer to different properties :-). How did you figure out that there is 0.14 lambda peak-to-peak OPD. Shouldn't this be 0.14 lambda peak-to-valley OPD

[jordens]

GeometricTrace(...).opd() will give you an OPD map. .psf() will give you the PSF map. These methods should be correct but they are not tested as thoroughly as some other parts. You get about 0.050 lambda RMS OPD in the paraxial focus and 0.035 in minimum geometric RMS focus.

s.paraxial.refocus()
g = ro.GeometricTrace(s)
g.rays_point((0, 0), nrays=20, distribution="radau", filter=False, clip=False)
# g.refocus()
p, q, opd = g.opd(resample=False)
np.sqrt((opd**2 * g.w).sum()/g.w.sum())

Peak-to-peak and peak-to-vally are the same to me (https://en.wikipedia.org/wiki/Amplitude#Peak-to-peak_amplitude)

[hstarmans]

Many thanks!! Really helped!! Please note that the model and Rayopt agree :+1:

lambda RMS optical path difference:
                               Rayopt                model 
10 degrees                0.006                  0.005
24 degrees                0.034                  0.032
30 degrees                0.054                  0.05