elena/gaussianoptics
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
gaussianoptics/gaussianoptics.py
Elena Arenskötter 47ea699f3e Add repo
2023-11-04 13:44:26 +01:00

255 lines
12 KiB
Python
Raw Blame History

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]
Reference in New Issue View Git Blame Copy Permalink