import numpy as np def gaussianoptics(opticalbench,lam,q0,startposition,N): ''' 'path', refractive index , position[mm], (relative: 'r', absolute: 'a' (position is the end position), filling: 'f') 'thinlens', f [mm] , position[mm], (relative: 'r', absolute: 'a') 'thicklens', [ n1 = refractive index outside of the lens, n2 = refractive index of the lens itself (inside the lens), R1 = Radius of curvature of First surface, R2 = Radius of curvature of Second surface, t = center thickness of lens] , position[mm], (relative: 'r', absolute: 'a') 'flatmirror', 0 , position[mm], (relative: 'r', absolute: 'a') 'curvedmirror', [R [mm], angle [°], plane (0=tangential, 1=sagittal)], position[mm], (relative: 'r', absolute: 'a') 'flatrefraction', [ n1 = initial refractive index, n2 = final refractive index] , position[mm], (relative: 'r', absolute: 'a') 'curvedrefraction', [ n1 = initial refractive index, n2 = final refractive index, R = radius of curvature R > 0 for convex (centre of curvature after interface)], position[mm], (relative: 'r', absolute: 'a') example: opticalbench = ['path',1,0.050,'a''thinlens',0.05,0.100,'r''path',[1,0,0,0,0],0,'f''thinlens',[0.15,0,0,0,0],0.200,'r'] TODO: - check if 'path' is first and last element - find 'empty space' - normalize curve by refractive index other than 1 Elena Arenskötter 05.08.2020 ''' elements = len(opticalbench) OB = [[0,0,0,0] for j in range(elements)] #[optical element, property, start position, stop position] position = startposition #find absolute positions for i in range(0,elements): if opticalbench[i][0] == 'path': OB[i][0] = 'path' OB[i][1] = opticalbench[i][1] #properties #start position OB[i][2] = position #stop position if opticalbench[i][3] == 'r': position = position + opticalbench[i][2] OB[i][3] = position elif opticalbench[i][3] == 'a': if opticalbench[i][2] < position: disp('Check sequential arrangement!') break else: position = opticalbench[i][2] OB[i][3] = position elif opticalbench[i][3] == 'f': if i == elements-1: position = position + 100 OB[i][3] = position print('Warning: Floating path at the end of your optical bench!') print('-> Replaced by 100mm relative length.') print('Optical bench entry: '+str(i)) else: if opticalbench[i+1][3] == 'a': position = opticalbench[i+1][2] OB[i][3] = position elif opticalbench[i+1][3] == 'r': position = position + opticalbench[i+1][2] OB[i][3] = position elif opticalbench[i+1][3] == 'f': print('Warning: Consecutive floating/relative elements!') print('Optical bench entry: '+str(i)) break else: print('Warning: Unknown position parameter.') print('Optical bench entry: ',str(i+1)) break else: print('Warning: Unknown position parameter.') print('Optical bench entry: '+str(i)) break else: if opticalbench[i][0] == 'thinlens': OB[i][0] = 'thinlens' elif opticalbench[i][0] == 'thicklens': OB[i][0] = 'thicklens' elif opticalbench[i][0] == 'flatmirror': OB[i][0] = 'flatmirror' elif opticalbench[i][0] == 'curvedmirror': OB[i][0] = 'curvedmirror' elif opticalbench[i][0] == 'flatrefraction': OB[i][0] = 'flatrefraction' elif opticalbench[i][0] == 'curvedrefraction': OB[i][0] = 'curvedrefraction' else: print('Warning: Unknown optical element!') print('--> '+opticalbench[i][0]+' entry: '+str(i)) break OB[i][1]=opticalbench[i][1] #properties #start position OB[i][2] = position #stop position if opticalbench[i][3] == 'a': if opticalbench[i][2] < position: print('Check sequential arrangement!') print('Optical bench entry: '+str(i)) break else: position = opticalbench[i][2] OB[i][3] = opticalbench[i][2] elif opticalbench[i][3] == 'r': if i == 0: position = position + opticalbench[i][2] OB[i][3] = position else: if opticalbench[i-1][3] == 'a': position = position + opticalbench[i][2] OB[i][3] = position elif opticalbench[i-1][3] == 'r': position = position + opticalbench[i][2] OB[i][3] = position elif opticalbench[i-1][3] == 'f': OB[i][3] = position else: print('Warning: Unknown position parameter.') print('Optical bench entry: '+str(i-1)) break else: print('Warning: Unknown position parameter.') print('Optical bench entry: '+str(i)) break #create ray-matrix M = [[0,0,0] for j in range(elements)] M_tmp = np.empty((2,2)) M_tmp[:] = 0 #[start position, stop position, matrix] parameter is always z for i in range(0,elements): M[i][0] = OB[i][2] M[i][1] = OB[i][3] if OB[i][0] == 'path': if i == 0: M[i][2] = 'np.array([[1,(z-'+str(OB[i][2])+')/'+str(OB[i][1])+'],[0,1]])' else: z = OB[i][2] M_tmp = eval(M[i-1][2]) M[i][2] = 'np.matmul(np.array([[1,(z-'+str(OB[i][2])+')/'+str(OB[i][1])+'],[0,1]]),\ np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']]))' elif OB[i][0] == 'thinlens': if i == 0: M[i][2] = 'np.array([[1,0],[-1/'+str(OB[i][1])+',1]])' else: z = OB[i][2] M_tmp = eval('np.matmul(np.array([[1,0],[-1/'+str(OB[i][1])+',1]]),'+M[i-1][2]+')') M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' elif OB[i][0] == 'thicklens': # [1,0;n2-n1/(R2*n1),n2/n1] * [1,d/n20,1] * [1,0;(n1-n2)/(R1*n2),n1/n2] M_tmp = np.matmul(np.array([[1,0],[(OB[i][1][1]-OB[i][1][0])/(OB[i][1][3]*OB[i][1][0]),OB[i][1][1]/OB[i][1][0]]]), np.matmul(np.array([[1,OB[i][1][4]/OB[i][1][1]],[0,1]]), np.array([[1,0],[(OB[i][1][0]-OB[i][1][1])/(OB[i][1][2]*OB[i][1][1]),OB[i][1][0]/OB[i][1][1]]]))) if i == 1: M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' else: z = OB[i][2] M_tmp = eval('np.matmul(np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']]),'+M[i-1][2]+')') M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' elif OB[i][0] == 'flatmirror': M_tmp = np.array([[1,0],[0,1]]) if i == 1: M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' else: z = OB[i][2] M_tmp = eval('np.matmul(np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']]),'+M[i-1][2]+')') M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' elif OB[i][0] == 'curvedmirror': if OB[i][1][2] == 0: #0=tangential, 1=sagittal M_tmp = np.array([[1,0],[-2/(OB[i][1][0]*np.cos(np.pi*OB[i][1][1]/180)),1]]) else: M_tmp = np.array([[1,0],[-2*np.cos(np.pi*OB[i][1][1]/180)/OB[i][1][0],1]]) if i == 1: M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' else: z = OB[i][2] M_tmp = eval('np.matmul(np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']]),'+M[i-1][2]+')') M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' elif OB[i][0] == 'flatrefraction': #M_tmp = [1,00,OB[i,2](2)/OB[i,2](1)] M_tmp = np.array([[1,0],[0,OB[i][1][0]/OB[i][1][1]]]) if i == 1: M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' else: z = OB[i][2] M_tmp = eval('np.matmul(np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']]),'+M[i-1][2]+')') M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' elif OB[i][0] == 'curvedrefraction': # [1 , 0 ] # [(n1-n2)/R*n2 , n1/n2] M_tmp = np.array([[1,0],[(OB[i][1][0]-OB[i][1][1])/(OB[i][1][2]*OB[i][1][1]),OB[i][1][0]/OB[i][1][1]]]) if i == 1: M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' else: z = OB[i][2] M_tmp = eval('np.matmul(np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']]),'+M[i-1][2]+')') M[i][2] = 'np.array([['+str(M_tmp[0][0])+','+str(M_tmp[0][1])+'],['+str(M_tmp[1][0])+','+str(M_tmp[1][1])+']])' x = np.linspace(M[0][0],M[-1][1],N) R = np.zeros((N,2)) w = np.zeros((N,1)) j = 0 for i in range(np.size(x)): z = x[i] while M[j][1] < z: j = j+1 M_tmp = eval(M[j][2]) q = ( M_tmp[0][0] * q0 + M_tmp[0][1] ) / ( M_tmp[1][0] * q0 + M_tmp[1][1]) if np.real(q) == 0: R[i] = 0 else: R[i] = np.real( q ) * (1 + ( np.imag( q ) / np.real( q ))**2) if OB[j][0] == 'path': w[i] = np.sqrt(lam*np.abs(np.imag(q))/np.pi)*np.sqrt(1+(np.real(q)/np.imag(q))**2)/np.sqrt(OB[j][1]) else: w[i] = np.sqrt(lam*np.abs(np.imag(q))/np.pi)*np.sqrt(1+(np.real(q)/np.imag(q))**2) M_tot = M_tmp return [x,w,R,OB,M_tot]