incremented version and minor updates

Author: ptbrown1729

scripts/2024_08_28_perfect_lens_pupil_translation.py

@@ -0,0 +1,124 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import raytrace.raytrace as rt
+from raytrace.materials import Vacuum, Ebaf11, Nsf11
+
+na = 0.98
+f = 1
+wlen = 0.561
+
+wlen = 0.635
+
+l1 = rt.Doublet(Ebaf11(),
+                Nsf11(),
+                radius_crown=50.8,
+                radius_flint=-247.7,
+                radius_interface=-41.7,
+                thickness_crown=20.,
+                thickness_flint=3.,
+                aperture_radius=25.4,
+                input_collimated=False,
+                names="AC508-075-A-ML")
+
+l2 = rt.Doublet(Ebaf11(),
+                Nsf11(),
+                radius_crown=50.8,
+                radius_flint=-247.7,
+                radius_interface=-41.7,
+                thickness_crown=20.,
+                thickness_flint=3.,
+                aperture_radius=25.4,
+                input_collimated=True,
+                names="AC508-075-A-ML")
+
+# flat surface at focal plane of lens 1
+cp1 = l1.get_cardinal_points(wlen, Vacuum(), Vacuum())
+f1_left = cp1[0][-1]
+f1_right = cp1[1][-1]
+wd_right = f1_right - l1.surfaces[-1].paraxial_center[-1]
+
+system = rt.System([rt.FlatSurface([0, 0, 0], [0, 0, 1], 25.4)],
+                   []
+                   )
+system = system.concatenate(l1, Vacuum(), -f1_left)
+
+# find collimated distance between l1 and l2
+d = l2.find_paraxial_collimated_distance(l2, wlen, Vacuum(), Vacuum(), Vacuum())
+
+# add flat surface at Fourier plane
+system = system.concatenate(rt.FlatSurface([0, 0, 0],
+                                           [0, 0, 1],
+                                           25.4),
+                            Vacuum(),
+                            wd_right)
+ind_pupil = len(system.surfaces) - 1
+
+# add lens #2
+system = system.concatenate(l2, Vacuum(), d - wd_right)
+
+
+# last surface at focal plane
+c2 = l2.get_cardinal_points(wlen, Vacuum(), Vacuum())
+wd2 = c2[1][2] - l2.surfaces[-1].paraxial_center[2]
+
+system = system.concatenate(rt.FlatSurface([0, 0, 0],
+                                           [0, 0, 1],
+                                           25.4),
+                            Vacuum(),
+                            wd2)
+
+# add lenses
+system_detect = rt.System([rt.FlatSurface([0, 0, 0],
+                                   [0, 0, 1],
+                                  f * na),
+                    rt.PerfectLens(f,
+                                  [0, 0, f],
+                                  [0, 0, 1],
+                                  np.arcsin(na)),
+                    rt.FlatSurface([0, 0, 2*f],
+                                   [0, 0, 1],
+                                   f * na
+                                   ),
+                    rt.PerfectLens(f,
+                                  [0, 0, 3*f],
+                                  [0, 0, 1],
+                                  np.arcsin(na)),
+                    rt.FlatSurface([0, 0, 4*f],
+                                   [0, 0, 1],
+                                   f * na
+                                   )
+                                  ],
+                    [Vacuum(), Vacuum(), Vacuum(), Vacuum()]
+                                  )
+
+system = system.concatenate(system_detect,
+                            Vacuum(),
+                            0
+                            )
+
+nrays = 101
+rays = rt.get_ray_fan([0, 0, 0],
+                   10 * np.pi/180,
+                   nrays,
+                   wlen)
+rays_out = system.ray_trace(rays, Vacuum(), Vacuum())
+
+system.plot(rays_out)
+plt.show()
+
+figh = plt.figure()
+grid = figh.add_gridspec(nrows=1, ncols=1)
+
+ax1 = figh.add_subplot(grid[0])
+ax1.set_title('pupil after')
+ax1.plot(rays_out[-1, :, 0],
+         rays_out[-1, :, 6] - rays_out[-1, nrays // 2, 6],
+         label="pupil after")
+ax1.plot(rays_out[-9, :, 0],
+         rays_out[-9, :, 6] - rays_out[-9, nrays //2 , 6],
+         'r.',
+         label='pupil before')
+ax1.set_xlabel("Height (mm)")
+ax1.set_ylabel("OPL (mm)")
+ax1.legend()
+plt.show()

src/raytrace/__init__.py

@@ -1 +1 @@
-__version__ = "0.0.0"
+__version__ = "0.1.0"

src/raytrace/raytrace.py

@@ -14,12 +14,11 @@
 from typing import Optional, Union
 from collections.abc import Sequence
 from numpy.typing import NDArray
-import copy
+from copy import deepcopy
 import numpy as np
 from matplotlib.figure import Figure
 from matplotlib.axes._axes import Axes
 import matplotlib.pyplot as plt
-import warnings
 from raytrace.materials import Material, Vacuum
 
 try:
@@ -69,17 +68,17 @@ def get_ray_fan(pt: NDArray,
     if np.linalg.norm(center_ray) != 1:
         raise ValueError("center_ray must be a unit vector")
 
-    thetas = np.linspace(-theta_max, theta_max, n_thetas)
-    phis = np.arange(nphis) * 2*np.pi / nphis
-    rays = np.zeros((n_thetas * nphis, 8))
+    thetas = xp.linspace(-theta_max, theta_max, n_thetas)
+    phis = xp.arange(nphis) * 2*np.pi / nphis
+    rays = xp.zeros((n_thetas * nphis, 8))
 
-    tts, pps = np.meshgrid(thetas, phis)
+    tts, pps = xp.meshgrid(thetas, phis)
     tts = tts.ravel()
     pps = pps.ravel()
 
-    enx = np.cross(np.array([0, 1, 0]), center_ray)
-    enx = enx / np.linalg.norm(enx)
-    eny = np.cross(center_ray, enx)
+    enx = xp.cross(np.array([0, 1, 0]), center_ray)
+    enx = enx / xp.linalg.norm(enx)
+    eny = xp.cross(center_ray, enx)
 
     pt = np.array(pt).squeeze()
 
@@ -405,7 +404,7 @@ def reverse(self):
         flip direction of the optic we are considering (so typically rays now enter from the right)
         :return:
         """
-        surfaces_rev = [copy.deepcopy(self.surfaces[-ii]) for ii in range(1, len(self.surfaces) + 1)]
+        surfaces_rev = [deepcopy(self.surfaces[-ii]) for ii in range(1, len(self.surfaces) + 1)]
 
         for ii in range(len(self.surfaces)):
             surfaces_rev[ii].input_axis *= -1
@@ -436,13 +435,13 @@ def concatenate(self,
 
         # specify distance between surfaces as distances between the paraxial foci
         if isinstance(other, System):
-            new_surfaces = [copy.deepcopy(s) for s in other.surfaces]
+            new_surfaces = [deepcopy(s) for s in other.surfaces]
             new_materials = other.materials
             other_stop = other.aperture_stop
             new_surfaces_by_name = other.surfaces_by_name
             new_names = other.names
         elif isinstance(other, Surface):
-            new_surfaces = [copy.deepcopy(other)]
+            new_surfaces = [deepcopy(other)]
             new_materials = []
             other_stop = None
             new_surfaces_by_name = np.array([0])
@@ -1284,9 +1283,9 @@ def propagate(self,
         ds_out = mag_na * normals + mag_nc * nc
 
         rays_refracted = xp.concatenate((rays_intersection[:, :3],
-                                   ds_out,
-                                   rays_intersection[:, 6:]),
-                                  axis=1)
+                                         ds_out,
+                                         rays_intersection[:, 6:]),
+                                         axis=1)
         rays_refracted[xp.isnan(ds_out[:, 0]), :3] = xp.nan
 
         # check array was within aperture