[Demo only, do not merge] Define a thin lens to study infinite conjugate systems

opticspy/ray_tracing/draw.py

@@ -14,19 +14,33 @@ def draw_surface(r,x0,d1,d2):
     verts_1 = []
     verts_2 = []
     codes = []
-    for y in np.linspace(0,d1/2,10):
-        if r > 0:
-            x = -np.sqrt(r**2-y**2) + r
-        else:
-            x = np.sqrt(r**2-y**2) + r
-        verts_1.append([x+x0,y])
-        verts_2.append([x+x0,-y])
+    r = (r if r is not None else np.inf)
+
+    # Sample an arc with 10 control points
+    ymax = (6.0 if np.isnan(d1) else d2) / 2
+    y = np.linspace(0, ymax, 10)
+
+    sign = (-1.0 if r > 0 else 1.0)
+    x = sign * np.sqrt(r**2 - y**2) + r
+
+    # Top half
+    verts_1 = np.vstack((x + x0, y)).T
+
+    # Bottom half
+    verts_2 = np.vstack((x + x0, -y)).T
+
     if d1 == d2:
-        verts = verts_1[::-1] + verts_2[1:]
+        # If the curved surface meets the lens edge, concat the vertices
+        verts = np.vstack((verts_1[::-1], verts_2[1:]))
     elif d1 < d2:
-        verts = [[verts_1[-1][0],d2/2]] + verts_1[::-1] + verts_2[1:] + [[verts_2[-1][0],-d2/2]]
+        # If the dome is smaller than the lens edge, add straight edges
+        verts = np.vstack(([verts_1[-1][0],d2/2], verts_1[::-1], verts_2[1:], [verts_2[-1][0],-d2/2]))
+    else:
+        print(f'Warning: ray position larger than lens diameter: d1: {d1:1.3g}, d2: {d2:1.3g}')
+        verts = np.vstack((verts_1[::-1], verts_2[1:]))
+
     codes.append(Path.MOVETO)
-    for j in range(len(verts)-1):
+    for j in range(verts.shape[0] - 1):
         codes.append(Path.LINETO)
     return verts,codes,[verts[0],verts[-1]]
 
@@ -117,6 +131,13 @@ def draw_system(Lens):
         patch = patches.PathPatch(path, facecolor='none', lw=1)
         ax.add_patch(patch)
 
+        # Also plot the control points
+        def draw_control_points(path):
+            x, y = zip(*path.vertices)
+            ax.plot(x, y, 'g+-')
+
+        #draw_control_points(path)
+
         # start drawing edge
         if num == 0:
             pass
@@ -146,8 +167,10 @@ def draw_system(Lens):
     ax.set_xlim(-0.5*sum(thinkness_list[1:-1]),1.1*sum(thinkness_list[1:-1]))
     d = max(draw_diameter_list)
     ax.set_ylim(-d/4*3,d/4*3)
+    ax.set_ylim(-1, 1)
     plt.axis('equal')
-    plt.show()
+    plt.xlabel('x / mm')
+    plt.ylabel('y / mm')
 
 
 def draw_rays(Lens):

opticspy/ray_tracing/first_order_tools.py

@@ -3,6 +3,7 @@
 import numpy as __np__
 import matplotlib.pyplot as __plt__
 
+FOCAL_LENGTH = 50.0
 
 def start_end(Lens,start_surface,end_surface):
 	'''
@@ -31,6 +32,9 @@ def R(c,n_left,n_right):
 	'''
 	return __np__.array([[1,0],[-c*(n_right-n_left),1]])
 
+def lens(focal_length, n_right):
+	return __np__.array([[1,0],[-1/focal_length,1/n_right]])
+
 def ABCD(matrix_list):
 	'''
 	ABCD matrix calculator
@@ -51,23 +55,28 @@ def ABCD_start_end(Lens,start_surface=0,end_surface=0):
 	'''
 	s = Lens.surface_list
 	start,end = start_end(Lens,start_surface,end_surface)
+
 	R_matrix = []
 	T_matrix = []
 	RT_matrix = []
-	for i in range(start,end+1):
-		i = i - 1
-		# print 'i:',i
-		# print 'surface_num',s[i].number
-		c = 1/s[i].radius
+	for i in range(start - 1,end):
 		# now use central wavelength as reference
 		n_left = s[i-1].indexlist[int(len(s[i-1].indexlist)/2)]
 		n_right = s[i].indexlist[int(len(s[i].indexlist)/2)]
-		t = s[i].thickness
+
+		if s[i].radius is None:
+			# Ideal thin lens
+			RT_matrix.append(lens(FOCAL_LENGTH,n_right))
+			assert s[i+1].thickness == FOCAL_LENGTH
+			continue
+
+		# print 'surface_num',s[i].number
+		c = 1/s[i].radius
 		R1 = R(c,n_left,n_right)
-		T1 = T(t,n_right)
-		# print R1
-		# print T1
 		RT_matrix.append(R1)
+
+		t = s[i].thickness
+		T1 = T(t,n_right)
 		if i+1 != end:
 			RT_matrix.append(T1)
 	# print '--------------------------------'

opticspy/ray_tracing/trace.py

@@ -3,6 +3,7 @@
 import matplotlib.pyplot as __plt__
 from . import field,first_order_tools,surface
 from . import output_tools
+from .first_order_tools import FOCAL_LENGTH
 
 # Function: trace rays
 # input a list of ray
@@ -231,8 +232,7 @@ def traceray(ray_list, surface1, surface2, wave_num):
         #print 'Pos',Pos
         KLM = ray.KLM
         #print 'KLM:',KLM
-        c1 = 1 / surface1.radius
-        c2 = 1 / surface2.radius
+        c2 = (None if surface2.radius is None else 1 / surface2.radius)
         n1 = surface1.indexlist[wave_num-1]
         n2 = surface2.indexlist[wave_num-1]
         tn1 = surface1.thickness
@@ -243,9 +243,17 @@ def traceray(ray_list, surface1, surface2, wave_num):
         Pos_new = xyz + delta * KLM
         Pos_new_list.append(Pos_new)
         # calculate new ray direction
-        # if curvature == 0, it is a stop, object or image plane
-        # don't need to calculate the new ray direction
-        if c2 == 0:
+        if c2 == None:
+            # Assume an ideal eye-piece having an effective focal length of 10mm
+            Kp = KLM[0] - Pos_new[0] / FOCAL_LENGTH
+            Lp = KLM[1] - Pos_new[1] / FOCAL_LENGTH
+            Mp = KLM[2] - Pos_new[2] / FOCAL_LENGTH
+            KLM_new = __np__.asarray([Kp, Lp, Mp])
+            #KLM_new = __np__.asarray(KLM)
+
+        elif c2 == 0:
+            # if curvature == 0, it is a stop, object or image plane
+            # don't need to calculate the new ray direction
             KLM_new = KLM
         else:
             sigma = __np__.sqrt(n2 ** 2 - n1 ** 2 * (1 - cosI ** 2)) - n1 * cosI
@@ -265,7 +273,7 @@ def pos(Pos, KLM, curvature):
     '''
     calculate new position of a ray on spherical surface
     '''
-    c = curvature
+    c = (0 if curvature is None else curvature)
     x0 = Pos[0]
     y0 = Pos[1]
     z0 = Pos[2]