diff --git a/gaussianoptics.py b/gaussianoptics.py index 25fe5c4..947511d 100644 --- a/gaussianoptics.py +++ b/gaussianoptics.py @@ -1,254 +1,254 @@ -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 - Jan 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] - - +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] + +