%pylab inline
import numpy as np
import matplotlib.pyplot as plt
from rayopt import *
np.set_printoptions(precision=3)
Populating the interactive namespace from numpy and matplotlib
description = "oslo cooke triplet example 50mm f/4 20deg"
columns = "type roc distance radius material"
text = """
O 0 0 .364 AIR
S 21.25 5 6.5 SK16
S -158.65 2 6.5 AIR
S -20.25 6 5 F4
S 19.3 1 5 AIR
A 0 0 4.75 AIR
S 141.25 6 6.5 SK16
S -17.285 2 6.5 AIR
I 0 42.95 .364 AIR
"""
_description = "triplet 50mm f/4 20deg"
_columns = "type curvature distance radius material"
_text = """
O 0 0 .364 AIR
S .25285 5 1.8 1.62
S -.01474 .6 1.8 AIR
S -.1994 1.0654 1.3 1.621
S .25973 .15 1.3 AIR
A 0 .1 1.1 AIR
S .05065 1.0396 1.7 1.62
S -.24588 .6 1.7 AIR
I 0 8.27937 2 AIR
"""
_description = "photo triplet, f/2.7, f=100 U.S.-Pat 2,453,260 (1948-Pestrecov)"
_columns = "type distance roc diameter material"
_text = """
O 0 0 .25 AIR
S 20 40.94 20 1.617/55
S 8.74 0 20 AIR
S 11.05 -55.65 20 1.649/33.8
S 2.78 39.75 20 AIR
A 0 0 30 AIR
S 7.63 107.56 30 1.617/55
S 9.54 -43.33 30 AIR
I 79.34 0 30 AIR
"""
_description="cooke type triplet, USP 2453260 Pestrecov, Modern Optical Engineering, Smith"
_columns = "type roc distance radius material"
_text = """
O 0 0 .3 AIR
S 40.94 10 16 S-BSM9
S 0 8.74 16 AIR
S -55.65 10.05 14 S-TIM22
A 0 2.78 12 S-TIM22
S 39.75 0 14 AIR
S 107.56 7.63 14.5 S-BSM9
S -43.33 9.54 14.5 AIR
I 0 90 0 AIR
"""
s = system_from_text(text, columns.split(),
description=description)
s.object.angle = np.deg2rad(20)
s = system_from_yaml("""
description: "code v cooke triplet example 50mm f/4.5 20deg"
stop: 3
object: {angle_deg: 20, pupil: {radius: 5.55}}
elements:
- {material: air, radius: 10}
- {distance: 4, material: SCHOTT/SK16, roc: 21.48138, radius: 7}
- {distance: 2, material: air, roc: -124.1, radius: 7}
- {distance: 5.26, material: HOYA/F4, roc: -19.1, radius: 4.1}
- {distance: 1.25, material: air, roc: 22, radius: 4.1}
- {distance: 4.69, material: SCHOTT/SK16, roc: 328.9, radius: 7}
- {distance: 2.25, material: air, roc: -16.7, radius: 7}
- {distance: 43.050484, material: air}
""")
s = system_from_yaml("""
description: "smith triplet p444, 100mm f/8 23.4deg"
stop: 4
object: {angle_deg: 23.4, pupil: {radius: 6.25}}
elements:
- {material: air, radius: 20}
- {distance: 10, material: SCHOTT/SK4, roc: 40.1, radius: 17}
- {distance: 6, material: air, roc: -537.1, radius: 17}
- {distance: 10, material: SCHOTT/FN11, roc: -47.0, radius: 15}
- {distance: 1, material: air, roc: 40, radius: 15}
- {distance: 10.8, material: SCHOTT/SK4, roc: 234.5, radius: 16}
- {distance: 6, material: air, roc: -37.9, radius: 16}
- {distance: 85.3, material: air, radius: 45}
""")
s.update()
s.paraxial.focal_length_solve(100)
s.paraxial.update()
s.paraxial.refocus()
print(s)
#s.reverse()
#print(s)
System: smith triplet p444, 100mm f/8 23.4deg Scale: 1.0 mm Wavelengths: 588, 656, 486 nm Fields: 0, 0.7, 1 Object: Semi-Angle: 23.4 deg Pupil: Pupil Distance: 28.5629 Refractive Index: 1.00028 Radius: 6.25 Image: Radius: 45 Update Radius: True Pupil: Pupil Distance: -104.811 Refractive Index: 1.00028 Update Radius: True Radius: 20.3004 Stop: 4 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 10 40.1 34 SCHOTT/SK4 1.613 1.613 58.63 2 S 6 -537.1 34 basic/air 1.000 1.000 89.30 3 S 10 -47 30 SCHOTT/FN11 1.621 1.621 36.18 4 S 1 40 30 basic/air 1.000 1.000 89.30 5 S 10.8 234.5 32 SCHOTT/SK4 1.613 1.613 58.63 6 S 6 -37.94 32 basic/air 1.000 1.000 89.30 7 S 85.48 inf 90 basic/air 1.000 1.000 89.30
a = Analysis(s, refocus_full=False)
/home/rj/work/nist/pyrayopt/rayopt/paraxial_trace.py:186: RuntimeWarning: divide by zero encountered in true_divide return self.n[(0, -2), ]/(2*na) /home/rj/work/nist/pyrayopt/rayopt/paraxial_trace.py:191: RuntimeWarning: divide by zero encountered in true_divide return 1.22*self.wavelength/(2*na)/self.system.scale
System: smith triplet p444, 100mm f/8 23.4deg Scale: 1.0 mm Wavelengths: 588, 656, 486 nm Fields: 0, 0.7, 1 Object: Semi-Angle: 23.4 deg Pupil: Pupil Distance: 28.5629 Refractive Index: 1.00028 Radius: 6.25 Image: Radius: 45 Update Radius: True Pupil: Pupil Distance: -104.984 Refractive Index: 1.00028 Update Radius: True Radius: 20.2931 Stop: 4 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 10 40.1 34 SCHOTT/SK4 1.613 1.613 58.63 2 S 6 -537.1 34 basic/air 1.000 1.000 89.30 3 S 10 -47 30 SCHOTT/FN11 1.621 1.621 36.18 4 S 1 40 30 basic/air 1.000 1.000 89.30 5 S 10.8 234.5 32 SCHOTT/SK4 1.613 1.613 58.63 6 S 6 -37.94 32 basic/air 1.000 1.000 89.30 7 S 85.48 inf 90 basic/air 1.000 1.000 89.30 lagrange: -2.7972 track length: 33.8 object, image height: [ 12.78 44.742] front, back focal length (from PP): [-100. 100.] entry, exit pupil height: [ 6.25 6.561] entry, exit pupil distance: [ 18.563 -19.504] front, back focal distance: [-76.668 85.457] front, back principal distance: [ 23.332 -14.543] front, back nodal distance: [ 23.332 -14.543] front, back numerical aperture: [ 0. 0.062] front, back f number: [ 8. 7.62] front, back working f number: [ inf 8.016] front, back airy radius: [ inf 0.006] transverse, angular magnification: [-0. 0.953] # T path n axial y axial nu chief y chief nu 0 S 0 1 6.25 0 -12.78 0.4475 1 S 10 1.613 6.25 -0.09546 -8.305 0.5744 2 S 16 1 5.895 -0.1022 -6.168 0.5814 3 S 26 1.621 4.873 -0.03782 -0.3558 0.5767 4 S 27 1 4.85 0.03744-2.165e-15 0.5767 5 S 37.8 1.613 5.254 0.02372 6.227 0.5605 6 S 43.8 1 5.343 -0.06252 8.312 0.4263 7 S 129.3 1 0 -0.06252 44.74 0.4263 # T SA3 CMA3 AST3 PTZ3 DIS3 TACHC TCHC 0 S 0 0 0 0 0 0 0 1 S -0.04459 -0.06875 -0.106 -0.5924 -1.077 -0.101 -0.1557 2 S -0.03949 0.2069 -1.084 -0.04423 5.913 -0.06911 0.3621 3 S 0.145 -0.4147 1.186 0.5097 -4.851 0.1699 -0.486 4 S 0.05063 0.184 0.6685 0.5989 4.605 0.1303 0.4736 5 S -0.00426 -0.04294 -0.4328 -0.1013 -5.384 -0.03258 -0.3284 6 S -0.1265 0.1289 -0.1313 -0.6261 0.7716 -0.1126 0.1147 7 S 0 -0 0 0 0 -0 0 -0.01928 -0.006555 0.1003 -0.2556 -0.02256 -0.01501 -0.01975
class GeomOp(Operand):
def __init__(self, system, rays=13, *args, **kwargs):
super(GeomOp, self).__init__(system, *args, **kwargs)
self.t = [GeometricTrace(self.system) for f in self.system.fields]
self.rays = rays
def get(self):
for t, f in zip(self.t, self.system.fields):
t.rays_point((0, f), nrays=self.rays,
distribution="radau",
clip=False, filter=False)
# for t, f in zip(self.t, self.system.fields):
# t.propagate()
v = np.concatenate([
(t.y[-1, 1:, :2] - t.y[-1, t.ref, :2])/np.sqrt(t.y.shape[1])
for t in self.t])
return np.where(np.isfinite(v), v, 10).ravel()
class PolyOp(Operand):
def __init__(self, system, kmax=3, *args, **kwargs):
super(PolyOp, self).__init__(system, *args, **kwargs)
self.t = PolyTrace(self.system, kmax=kmax)
def get(self):
self.t.propagate()
return self.t.stvwof[-1, :2, 1:].ravel()
s1 = system_from_yaml(system_to_yaml(s))
s1[1].material = s[5].material = Library.one().get("material", "LAK33")
s1[3].material = Library.one().get("material", "BASF2")
#s1.stop = 3
#s1.aperture.radius *= 1.5
s1.object.pupil.radius *= 1.5
s1.update()
variables = []
variables.extend(PathVariable(s1, (i, "distance"), (1, 6))
for i in (2, 4, 6))
variables.extend(PathVariable(s1, (i, "distance"), (.1, 15))
for i in (3, 5))
variables.extend(PathVariable(s1, (i, "curvature"), (-1/20, 1/20))
for i in (1, 2, 3, 4, 5))
def get(s):
s.update()
s.paraxial.focal_length_solve(100)
s.paraxial.propagate(start=-2)
s.paraxial.refocus()
s.paraxial.propagate(start=-1)
s.paraxial.resize()
s.paraxial.aberrations()
operands = [
FuncOp(s1, get),
FuncOp(s1, lambda s: s[-1].offset[2:], min=60),
FuncOp(s1, lambda s: s.edge_thickness()[1:], min=2),
#FuncOp(s1, lambda s: np.diff(s.track), min=2),
FuncOp(s1, lambda s: s.paraxial.transverse3[:, 5:].sum(0)/.5, min=-1, max=1),
FuncOp(s1, lambda s: s.paraxial.transverse3[:, 3:5].sum(0)/1, min=-1, max=1),
FuncOp(s1, lambda s: s.paraxial.transverse3[:, :3].sum(0)/.03, min=-1, max=1),
#FuncOp(s1, lambda s: s.paraxial.transverse3[:, (0, 1, 2, 5, 6)].sum(0), weight=1),
GeomOp(s1, rays=13, weight=1),
#PolyOp(s1, weight=1),
]
r = optimize(variables, operands, trace=True, tol=1e-4, options=dict( \
maxiter=100, disp=True, eps=1e-5, iprint=2))
r.accept()
print(s1)
NIT FC OBJFUN GNORM 1 12 1.034959E+00 1.291763E+02 2 25 4.422709E-01 3.357569E+01 3 38 3.455998E-01 2.530575E+01 4 51 1.119933E-01 1.402083E+01 5 64 2.167215E-02 3.952190E+00 6 77 6.887718E-03 1.222249E+00 7 90 6.169202E-03 1.017091E+00 8 103 5.508261E-03 1.035621E+00 9 115 4.507479E-03 4.204897E-01 10 128 4.378994E-03 4.886541E-01 11 140 4.259370E-03 3.619366E-01 12 152 4.112810E-03 3.850290E-01 13 164 3.669951E-03 3.262682E-01 14 176 2.959347E-03 1.970480E-01 15 188 3.119201E-03 6.546495E-01 16 200 3.012921E-03 4.186709E-01 17 212 2.969960E-03 3.887363E-01 Optimization terminated successfully. (Exit mode 0) Current function value: 0.00294449881868 Iterations: 17 Function evaluations: 213 Gradient evaluations: 17 System: smith triplet p444, 100mm f/8 23.4deg Scale: 1.0 mm Wavelengths: 588, 656, 486 nm Fields: 0, 0.7, 1 Object: Semi-Angle: 23.4 deg Pupil: Pupil Distance: 25.5827 Refractive Index: 1.00028 Radius: 9.375 Image: Radius: 45 Pupil: Pupil Distance: -93.5084 Refractive Index: 1.00028 Update Radius: True Radius: 8.76199 Stop: 4 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 10 32.88 32.694 SCHOTT/LAK33 1.754 1.754 52.43 2 S 5.9373 184.6 27.127 basic/air 1.000 1.000 89.30 3 S 7.7749 -82.93 15.205 SCHOTT/BASF2 1.664 1.664 35.83 4 S 1 30.96 14.357 basic/air 1.000 1.000 89.30 5 S 8.7232 110.2 25.118 SCHOTT/SK4 1.613 1.613 58.63 6 S 4.6083 -53.7 28.245 basic/air 1.000 1.000 89.30 7 S 79.321 inf 89.484 basic/air 1.000 1.000 89.30
fig, ax = plt.subplots(1, 3, figsize=(15, 7))
ax[0].plot(r.trace_x)
for i in range(len(r.trace_v[0])):
ax[1].plot([_[i] for _ in r.trace_v])
for i, _ in r.trace_f:
ax[2].semilogy(np.square(_).sum(-1), label=str(operands[i]))
ax[2].legend()
<matplotlib.legend.Legend at 0x7f7609d21668>
Analysis(s1, refocus_full=False)
/home/rj/work/nist/pyrayopt/rayopt/paraxial_trace.py:186: RuntimeWarning: divide by zero encountered in true_divide return self.n[(0, -2), ]/(2*na) /home/rj/work/nist/pyrayopt/rayopt/paraxial_trace.py:191: RuntimeWarning: divide by zero encountered in true_divide return 1.22*self.wavelength/(2*na)/self.system.scale /home/rj/work/nist/pyrayopt/rayopt/elements.py:208: RuntimeWarning: invalid value encountered in less_equal good = np.square(y[:, :2]).sum(1) <= self.radius**2
System: smith triplet p444, 100mm f/8 23.4deg Scale: 1.0 mm Wavelengths: 588, 656, 486 nm Fields: 0, 0.7, 1 Object: Semi-Angle: 23.4 deg Pupil: Pupil Distance: 25.5827 Refractive Index: 1.00028 Radius: 9.375 Image: Radius: 45 Pupil: Pupil Distance: -93.2576 Refractive Index: 1.00028 Update Radius: True Radius: 8.7429 Stop: 4 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 10 32.88 32.694 SCHOTT/LAK33 1.754 1.754 52.43 2 S 5.9373 184.6 27.127 basic/air 1.000 1.000 89.30 3 S 7.7749 -82.93 15.205 SCHOTT/BASF2 1.664 1.664 35.83 4 S 1 30.96 14.357 basic/air 1.000 1.000 89.30 5 S 8.7232 110.2 25.118 SCHOTT/SK4 1.613 1.613 58.63 6 S 4.6083 -53.7 28.245 basic/air 1.000 1.000 89.30 7 S 79.321 inf 89.484 basic/air 1.000 1.000 89.30 lagrange: -4.1957 track length: 28.044 object, image height: [ 11.446 44.742] front, back focal length (from PP): [-100. 100.] entry, exit pupil height: [ 9.375 8.743] entry, exit pupil distance: [ 15.583 -13.937] front, back focal distance: [-91.622 79.299] front, back principal distance: [ 8.378 -20.701] front, back nodal distance: [ 8.378 -20.701] front, back numerical aperture: [ 0. 0.093] front, back f number: [ 5.333 5.719] front, back working f number: [ inf 5.357] front, back airy radius: [ inf 0.004] transverse, angular magnification: [-0. 1.072] # T path n axial y axial nu chief y chief nu 0 S 0 1 9.375 0 -11.45 0.4475 1 S 10 1.754 9.375 -0.2149 -6.972 0.6074 2 S 15.94 1 8.648 -0.1796 -4.916 0.5873 3 S 23.71 1.664 7.252 -0.1215 -0.3512 0.5845 4 S 24.71 1 7.179 0.0325 2.387e-15 0.5845 5 S 33.44 1.613 7.462 -0.008958 5.097 0.5562 6 S 38.04 1 7.436 -0.09378 6.686 0.4799 7 S 117.4 1 0 -0.09378 44.74 0.4799 # T SA3 CMA3 AST3 PTZ3 DIS3 TACHC TCHC 0 S 0 0 0 0 0 0 0 1 S -0.284 -0.2345 -0.1936 -1.226 -1.172 -0.2338 -0.193 2 S -0.08908 0.3763 -1.589 0.2184 5.789 -0.1003 0.4237 3 S 0.374 -0.8285 1.835 0.4515 -5.065 0.23 -0.5095 4 S 0.2043 0.4516 0.9982 1.21 4.88 0.2255 0.4984 5 S -0.01435 -0.09032 -0.5686 -0.3232 -5.614 -0.05164 -0.3251 6 S -0.1932 0.2955 -0.452 -0.6637 1.707 -0.1193 0.1825 7 S 0 -0 0 0 0 -0 0 -0.002262 -0.02988 0.02999 -0.3338 0.5252 -0.04953 0.07705
<rayopt.analysis.Analysis at 0x7f7609a332e8>
from rayopt import GaussianTrace
s.object = FiniteConjugate(radius=.5)
g = GaussianTrace(s)
g.is_proper()
print(g.m, g.eigenvalues, g.real, g.stable)
g.refocus()
#qi = g.qi[0]
#g.qi[0, 0] *= 2
#qi[0, 1] = qi[1, 0] = 1e-2-1e-5j
#g.rays(qi)
#g.propagate()
#print g.qi[-1]
r = g.spot_radius_at(None)
print(r)
print(g.spot_radius_at(z=None, normal=True))
#print 1/g.qi.real
#print g.eigenmodes
#g.rays(np.eye(2)*g.eigenmodes[0])
g.propagate()
#print g.qi
print(g)
print(g.spot_radius)
fig, ax = plt.subplots()
s.plot(ax)
g.plot(ax, color="red", waist=True, scale=10)
#ax.set_xlim(0, 70)
#ax.set_ylim(-5, .5)
ax.set_aspect("auto")
plt.show()
[ 0.09 0.09] (array([ 0.09+0.996j, 0.09+0.996j]), array([ 0.09-0.996j, 0.09-0.996j])) [ True True] [ True True] [[ 0.5 0.5 ] [ 0.5 0.5 ] [ 0.472 0.472] [ 0.39 0.39 ] [ 0.388 0.388] [ 0.421 0.421] [ 0.426 0.426] [ 0.027 0.027]] (array([[ 0.5 , 0.5 ], [ 0.5 , 0.5 ], [ 0.472, 0.472], [ 0.39 , 0.39 ], [ 0.388, 0.388], [ 0.421, 0.421], [ 0.426, 0.426], [ 0.027, 0.027]]), array([ 0., 0., 0., 0., 0., 0., 0., 0.])) # T path spot a spot b spot ang waistx dz waisty dz waist x waist y 0 S 0 0.5 0.5 0 0 0 0.5 0.5 1 S 10 0.5 0.5 0 105.4 105.4 0.02447 0.02447 2 S 16 0.4716 0.4716 0 57.6 57.6 0.02286 0.02286 3 S 26 0.39 0.39 0 204.4 204.4 0.06124 0.06124 4 S 27 0.3881 0.3881 0 -125.9 -125.9 0.06143 0.06143 5 S 37.8 0.4206 0.4206 0 -414.6 -414.6 0.1088 0.1088 6 S 43.8 0.4263 0.4263 0 62.39 62.39 0.02742 0.02742 7 S 106.2 0.02742 0.02742 0-3.333e-14-1.218e-14 0.02742 0.02742 [[ 0.5 0.5 ] [ 0.5 0.5 ] [ 0.472 0.472] [ 0.39 0.39 ] [ 0.388 0.388] [ 0.421 0.421] [ 0.426 0.426] [ 0.027 0.027]]
/home/rj/work/nist/pyrayopt/rayopt/elements.py:570: RuntimeWarning: divide by zero encountered in true_divide s = .5*y*n0*(1 - mu)/l*(i + u/n) /home/rj/work/nist/pyrayopt/rayopt/elements.py:570: RuntimeWarning: invalid value encountered in true_divide s = .5*y*n0*(1 - mu)/l*(i + u/n) /home/rj/work/nist/pyrayopt/rayopt/elements.py:571: RuntimeWarning: invalid value encountered in double_scalars w = 4*k*n*(1 - mu)/l /home/rj/work/nist/pyrayopt/rayopt/elements.py:583: RuntimeWarning: invalid value encountered in double_scalars dc = s[1]*i[0]*i[1] + .5*(u[1]**2/n**2 - u0[1]**2/n0**2) + \ /home/rj/work/nist/pyrayopt/rayopt/elements.py:586: RuntimeWarning: invalid value encountered in true_divide tachc, tchc = -y[0]*i/l*(v0 - mu*v)