From 47ea699f3e2f42121e8615eb0ce4ff4e997ac3a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Elena=20Arensk=C3=B6tter?= Date: Sat, 4 Nov 2023 13:44:26 +0100 Subject: [PATCH] Add repo --- gaussianoptics.py | 254 ++++++++++++++++++++++++++++++++++++++++++++++ ionBeam854nm.py | 76 ++++++++++++++ test.py | 114 +++++++++++++++++++++ 3 files changed, 444 insertions(+) create mode 100644 gaussianoptics.py create mode 100644 ionBeam854nm.py create mode 100644 test.py diff --git a/gaussianoptics.py b/gaussianoptics.py new file mode 100644 index 0000000..25fe5c4 --- /dev/null +++ b/gaussianoptics.py @@ -0,0 +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] + + diff --git a/ionBeam854nm.py b/ionBeam854nm.py new file mode 100644 index 0000000..c5ef511 --- /dev/null +++ b/ionBeam854nm.py @@ -0,0 +1,76 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Aug 14 14:17:54 2020 + +@author: Jan +""" + + +import numpy as np +import matplotlib.pyplot as plt +import gaussianoptics as go + +#Laser wavelength +lam = 854e-9 + +#Number of points +N =5001 + +#Focal length and position definition +lens_coupler = 12e-3 +f1_tele = 30e-3 +pos_f1 = 100e-3 + +f2_tele = 300e-3 +pos_f2 = 330e-3 + +haloPosDif393_854 = 1.4e-3 + +f_halo1 = 25e-3 +pos_halo1 = 300e-3 + +f_halo2 = 25e-3 +pos_halo2 = 50e-3 - haloPosDif393_854 + +f1_colTele = 300e-3; +pos_colTele1= 400e-3; + +f2_colTele = 25e-3; +pos_colTele2= 497.5e-3; + +#Definition of optical arangement +opticalbench = [['path',1,0,'f'], + ['thinlens',12e-3,lens_coupler,'r'], + ['path',1,0,'f'], + ['thinlens',f1_tele,pos_f1,'r'], + ['path',1,0,'f'], + ['thinlens',f2_tele,pos_f2,'r'], + ['path',1,0,'f'], + ['thinlens',f_halo1,pos_halo1,'r'], + ['path',1,0,'f'], + ['thinlens',f_halo2,pos_halo2,'r'], + ['path',1,0,'f'], + ['thinlens',f1_colTele,pos_colTele1,'r'], + ['path',1,0,'f'], + ['thinlens',f2_colTele,pos_colTele2,'r'], + ['path',1,0.3,'r'] +] + + + +#Beamwaist inside fiber 780HP +w0 = 2.5e-6 +z0 = np.pi*w0**2/lam +q0 = complex(0,z0) + +[x,w,R,OB,M_tot] = go.gaussianoptics(opticalbench,lam,q0,0,N) + +fig, ax = plt.subplots() +ax.plot(x*1000, w*1000, color='black') +ax.plot(x*1000, -w*1000, color='black') + +plt.grid(True,which="both") + +plt.show() + + diff --git a/test.py b/test.py new file mode 100644 index 0000000..5844df6 --- /dev/null +++ b/test.py @@ -0,0 +1,114 @@ + +import numpy as np +import matplotlib.pyplot as plt +import gaussianoptics as go + + +lens_coupler = 12e-3 + + +Radius = 0.05 +alpha = 15 +R_mirr_t = -Radius*np.cos(alpha*np.pi/180) +R_mirr_s = -Radius/np.cos(alpha*np.pi/180) +pos_colTele2= 497.5e-3 + +opticalbench_tangential_trans = [['path',1,0,'f'], + ['thinlens',12e-3,lens_coupler,'r'], + ['path',1,0,'f'], + ['thicklens',[1,1.52,13.1e-3,np.inf,11.7e-3],20e-3,'r'], + ['path',1,0,'f'], + ['thinlens',50e-3,50e-3,'r'], + ['path',1,0,'f'], + ['flatmirror',0,50e-3,'r'], + ['path',1,0,'f'], + ['curvedmirror',[0.05,alpha,0],50e-3,'r'], + ['path',1,0,'f'], + ['flatrefraction',[1,1.52],50e-3,'r'], + ['path',1,0,'f'], + ['curvedrefraction',[1.52,1,R_mirr_t],15e-3,'r'], + ['path',1,0,'f'], + ['thinlens',100e-3,50e-3,'r'], + ['path',1,0.3,'r'] + ] + +opticalbench_sagital_trans = [['path',1,0,'f'], + ['thinlens',12e-3,lens_coupler,'r'], + ['path',1,0,'f'], + ['thicklens',[1,1.52,13.1e-3,np.inf,11.7e-3],20e-3,'r'], + ['path',1,0,'f'], + ['thinlens',50e-3,50e-3,'r'], + ['path',1,0,'f'], + ['flatmirror',0,50e-3,'r'], + ['path',1,0,'f'], + ['curvedmirror',[0.05,alpha,1],50e-3,'r'], + ['path',1,0,'f'], + ['flatrefraction',[1,1.52],50e-3,'r'], + ['path',1,0,'f'], + ['curvedrefraction',[1.52,1,R_mirr_s],15e-3,'r'], + ['path',1,0,'f'], + ['thinlens',100e-3,50e-3,'r'], + ['path',1,0.3,'r']] +''' +opticalbench = [['path',1,0,'f'], + ['thinlens',12e-3,lens_coupler,'r'], + ['path',1,0,'f'], + ['thinlens',f1_tele,pos_f1,'r'], + ['path',1,0,'f'], + ['thinlens',f2_tele,pos_f2,'r'], + ['path',1,0,'f'], + ['thinlens',f_halo1,pos_halo1,'r'], + ['path',1,0,'f'], + ['thinlens',f_halo2,pos_halo2,'r'], + ['path',1,0,'f'], + ['thinlens',f1_colTele,pos_colTele1,'r'], + ['path',1,0,'f'], + ['thinlens',f2_colTele,pos_colTele2,'r'], + ['path',1,0.3,'r'] +] +''' +lam = 854e-9 + +N =5001 + +#Beamwaist inside fiber 780HP +w0 = 2.5e-6 +z0 = np.pi*w0**2/lam +q0 = np.complex(0,z0) + + +fig, ax = plt.subplots() + +[x,w,R,OB,M_tot] = go.gaussianoptics(opticalbench_tangential_trans,lam,q0,0,N) +ax.plot(x*1000, w*1000, color='red') +ax.plot(x*1000, -w*1000, color='red') + +[x,w,R,OB,M_tot] = go.gaussianoptics(opticalbench_sagital_trans,lam,q0,0,N) +ax.plot(x*1000, w*1000, color='blue') +ax.plot(x*1000, -w*1000, color='blue') + +ax.set(xlabel='position (mm)', ylabel='beam radius (mm)', title='Testbench') + +plt.grid(True,which="both") + +Ymin = -10 +Ymax = 10 +elements = len(OB) +for i in range(elements): + if OB[i][0] == 'flatmirror': + ax.plot([OB[i][3]*1000, OB[i][3]*1000],[Ymin/1.1, Ymax/1.1],color='blue',linewidth=1) + elif OB[i][0] == 'path': + ax.plot([OB[i][3]*1000, OB[i][3]*1000],[Ymin/1.1, Ymax/1.1],color=[0.5,0.5,0.5],linewidth=1) + elif OB[i][0] == 'flatrefraction': + ax.plot([OB[i][3]*1000, OB[i][3]*1000],[Ymin/1.1, Ymax/1.1],color='red',linewidth=1) + elif OB[i][0] == 'thicklens': + ax.plot([OB[i][3]*1000, OB[i][3]*1000],[Ymin/1.1, Ymax/1.1],color='green',linewidth=1) + elif OB[i][0] == 'thinlens': + ax.plot([OB[i][3]*1000, OB[i][3]*1000],[Ymin/1.1, Ymax/1.1],color=[0.75,0,1],linewidth=1) + else: + ax.plot([OB[i][3]*1000, OB[i][3]*1000],[Ymin/1.1, Ymax/1.1],color='black',linewidth=1) + + +plt.show() + +