|
| 1 | +import numpy as np |
| 2 | +from LightPipes import * |
| 3 | +import matplotlib.pyplot as plt |
| 4 | +from matplotlib.offsetbox import OffsetImage, AnnotationBbox |
| 5 | +import matplotlib.image as mpimg |
| 6 | + |
| 7 | +''' |
| 8 | +Generation of an Airy beam. |
| 9 | +=========================== |
| 10 | +LightPipes for Python, Fred van Goor, 11-1-2022 |
| 11 | +
|
| 12 | +A non-diffracting Airy beam can be generated by substituting a cubic phase distribution |
| 13 | +in a laser beam. After a 2D Fourier transform with a positive lens the |
| 14 | +Airy beam will exist. |
| 15 | +With the parameters below, the Airy beam exists from z = 0 cm to z= 30 cm. |
| 16 | +Ref: |
| 17 | +G.A. Siviloglou, J. Broky, A. Dogariu and D.N. Christodoulides, PRL 99,213901(2007) |
| 18 | +T. Latychevskaia, D. Schachtler, and H.W. Fink, Applied Optics Vol. 55, Issue 22, pp. 6095-6101 (2016) |
| 19 | +''' |
| 20 | + |
| 21 | +wavelength = 650*nm #wavelength |
| 22 | +N=624 |
| 23 | +size=N*32.0*um |
| 24 | +dN1=N//2-40 |
| 25 | +dN2=N//2+40 |
| 26 | +beta=117/m # |
| 27 | +w0=10*mm #beam waist laser |
| 28 | +f=80*cm #focal length of Fourier lens |
| 29 | +k=2*np.pi/wavelength |
| 30 | + |
| 31 | +#Generate input beam |
| 32 | +F0=Begin(size,wavelength,N) |
| 33 | +F0=GaussHermite(F0,w0) |
| 34 | + |
| 35 | +#Cubic Phase Plate (CPP): |
| 36 | +X,Y = F0.mgrid_cartesian |
| 37 | +c = 2 * np.pi * beta |
| 38 | +c3 = (c**3)/3 |
| 39 | +CPP=c3*(X**3 + Y**3) + np.pi |
| 40 | +F0=SubPhase(F0,CPP) |
| 41 | +phase=Phase(F0) |
| 42 | + |
| 43 | +#initiate figure 1: |
| 44 | +fig1, axs1 = plt.subplots(nrows=3, ncols=1,figsize=(5.0,6.0)) |
| 45 | +figure1 = mpimg.imread('figure1.jpg') |
| 46 | +imagebox = OffsetImage(figure1, zoom=0.15) |
| 47 | +ab = AnnotationBbox(imagebox, (0.4, 0.6),frameon=False) |
| 48 | +s1=r'Generation of an Airy beam.'+'\n\n' |
| 49 | +s2=r'References:'+ '\n'\ |
| 50 | + r'(1) G.A. Siviloglou, J. Broky, A. Dogariu and D.N. Christodoulides, Phys. Rev. Lett, 99,213901(2007)'+'\n'\ |
| 51 | + r'(2) T. Latychevskaia, D. Schachtler, and H.W. Fink, Appl. Optics, 55, 6095-6101(2016)'+'\n\n' |
| 52 | +s3 = r'LightPipes for Python,' + '\n' + 'AiryBeam.py' + '\n\n'\ |
| 53 | + r'Parameters from ref 2:' + '\n'\ |
| 54 | + r'$\lambda = {:4.1f}$'.format(wavelength/nm) + r' $nm$' + '\n'\ |
| 55 | + r'$size = {:4.2f}$'.format(size/mm) + r' $mm$' + '\n'\ |
| 56 | + r'$N = {:4d}$'.format(N) + '\n'\ |
| 57 | + r'$w_0 = {:4.2f}$'.format(w0/mm) + r' $mm$'+ '\n'\ |
| 58 | + r'$f = {:4.2f}$'.format(f/cm) + r' $cm$'+ '\n'\ |
| 59 | + r'$\beta = {:4.2f}$'.format(beta/m) + r' $m^{-1}$'+ '\n'\ |
| 60 | + r'${\copyright}$ Fred van Goor, January 2022'+ '\n' |
| 61 | +axs1[0].text(0.,1.3,s1,fontsize=15, verticalalignment='top') |
| 62 | +axs1[0].text(0.,1.0,s2+s3,fontsize=5, verticalalignment='top');axs1[0].axis('off') |
| 63 | +axs1[1].add_artist(ab);axs1[1].axis('off') |
| 64 | +axs1[1].text(0.0,1.1,'figure 1 from reference 2',fontsize=5, verticalalignment='top') |
| 65 | +s=r'Phase distribution SLM' |
| 66 | +axs1[2].imshow(phase,cmap='gray'); axs1[2].axis('off'); axs1[2].set_title(s) |
| 67 | + |
| 68 | +#Fourier transform lens: |
| 69 | +F0=Lens(F0,f) |
| 70 | + |
| 71 | +#Propagate a distance f: |
| 72 | +F0=Fresnel(F0,f) |
| 73 | + |
| 74 | +def AiryBeam(z): |
| 75 | + ''' |
| 76 | + Propagates the Airy beam a distance z from the focal point.. |
| 77 | + returns a tuple: |
| 78 | + Position of the maximum intensity, |
| 79 | + The intensity distribution at a distance z, |
| 80 | + The peak intensity of the Airy beam at z. |
| 81 | + ''' |
| 82 | + F=Fresnel(F0,z) |
| 83 | + I=Intensity(F) |
| 84 | + #Find the coordinates of the maximum intensity: |
| 85 | + result = np.where(I == np.amax(I)) |
| 86 | + Coordinates = list(zip(result[0], result[1])) |
| 87 | + for Coord in Coordinates: |
| 88 | + Xmax=X[Coord[0],Coord[1]]; Ymax=Y[Coord[0],Coord[1]] |
| 89 | + print(Xmax/mm,Ymax/mm) #Xmax should be equal to Ymax |
| 90 | + return Xmax,Ymax,I,I[result][0] |
| 91 | + |
| 92 | +def xmax(z): |
| 93 | + ''' |
| 94 | + Returns: |
| 95 | + The theoretical deflection at z, |
| 96 | + The peak intensity of the Airy beam at a distance z, |
| 97 | + according to Latychevskaia et al. |
| 98 | + ''' |
| 99 | + c=z*f/(f+z) |
| 100 | + zdivc=(f+z)/f |
| 101 | + xmax=(1.0188*(beta*wavelength*f)-c**2/(4*k**2*(beta*wavelength*f)**3))*zdivc |
| 102 | + return xmax, (f/(f+z))**2 |
| 103 | + |
| 104 | +zstart=0*cm |
| 105 | +zend=31*cm |
| 106 | +z=np.arange(zstart,zend,1*cm) |
| 107 | +n=z.shape[0] |
| 108 | +deflection=np.zeros(n) |
| 109 | +Imax=np.zeros(n) |
| 110 | + |
| 111 | +#initiate figure 2: |
| 112 | +fig2, axs2 = plt.subplots(nrows=4, ncols=2,figsize=(6.0,6.0)) |
| 113 | +fig2.suptitle('simulation with data Latychevskaia et al.') |
| 114 | +fig2.subplots_adjust(hspace=0.8) |
| 115 | +for i in range(0,n): |
| 116 | + Xmax,Ymax,I,_Imax=AiryBeam(z[i]) |
| 117 | + |
| 118 | + if z[i] == 0*cm: |
| 119 | + axs2[0,0].imshow(I,cmap='jet') |
| 120 | + s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$' |
| 121 | + axs2[0,0].imshow(I,cmap='jet'); axs2[0,0].axis('off'); axs2[0,0].set_title(s) |
| 122 | + axs2[0,0].set_xlim(dN1,dN2) |
| 123 | + axs2[0,0].set_ylim(dN2,dN1) |
| 124 | + |
| 125 | + if z[i]== 10*cm: |
| 126 | + s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$' |
| 127 | + axs2[1,0].imshow(I,cmap='jet'); axs2[1,0].axis('off'); axs2[1,0].set_title(s) |
| 128 | + axs2[1,0].set_xlim(dN1,dN2) |
| 129 | + axs2[1,0].set_ylim(dN2,dN1) |
| 130 | + |
| 131 | + if z[i]== 20*cm: |
| 132 | + s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$' |
| 133 | + axs2[2,0].imshow(I,cmap='jet'); axs2[2,0].axis('off'); axs2[2,0].set_title(s) |
| 134 | + axs2[2,0].set_xlim(dN1,dN2) |
| 135 | + axs2[2,0].set_ylim(dN2,dN1) |
| 136 | + |
| 137 | + if z[i]== 30*cm: |
| 138 | + s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$' |
| 139 | + axs2[3,0].imshow(I,cmap='jet'); axs2[3,0].axis('off'); axs2[3,0].set_title(s) |
| 140 | + axs2[3,0].set_xlim(dN1,dN2) |
| 141 | + axs2[3,0].set_ylim(dN2,dN1) |
| 142 | + |
| 143 | + deflection[i]=Xmax |
| 144 | + Imax[i]=_Imax |
| 145 | + |
| 146 | +gs = axs2[0,1].get_gridspec() |
| 147 | +for ax in axs2[0:,-1]: |
| 148 | + ax.remove() |
| 149 | +ax1=fig2.add_subplot(gs[0:2,1]) |
| 150 | +ax1.plot(z/cm,-deflection/mm,'ro',z/cm,-xmax(z)[0]/mm) |
| 151 | +ax1.set_xlabel('z/cm'); ax1.set_ylabel('deflection/mm') |
| 152 | +ax1.legend(('simulated with LightPipes for Python', 'theory ref 2'), |
| 153 | + shadow=True, loc=(0.1, 0.8), handlelength=1.5, fontsize=5) |
| 154 | + |
| 155 | +ax2=fig2.add_subplot(gs[2:,1]) |
| 156 | +ax2.plot(z/cm,Imax/Imax[0],'ro',z/cm,xmax(z)[1]) |
| 157 | +ax2.set_xlabel('z/cm'); ax2.set_ylabel('Peak intensity') |
| 158 | + |
| 159 | +#Experimental data from ref 2, figure 2f: |
| 160 | +#In pixels, as good as possible.... |
| 161 | +X_Exp=np.array([122,151,184,213,245,276,308,336,366,397,429,458,491,520,551,582]) |
| 162 | +Y_Exp=np.array([106,139,177,216,198,207,247,266,289,327,336,354,350,360,360,366]) |
| 163 | +Imax_Exp=1.0-1.0/(428-106)*(Y_Exp-106) |
| 164 | +z_Exp=30*cm/(582-122)*(X_Exp-122) |
| 165 | +ax2.plot(z_Exp/cm,Imax_Exp,'+') |
| 166 | +ax2.legend(('simulated with LightPipes for Python', 'theory ref 2', 'experiment ref 2'), |
| 167 | + shadow=True, loc=(0.25, 0.8), handlelength=1.5, fontsize=5) |
| 168 | + |
| 169 | +plt.show() |
0 commit comments