From cc8238f64a036aad919ab004a55f66eebc1914ac Mon Sep 17 00:00:00 2001 From: salazarsilva Date: Thu, 24 Feb 2022 14:18:17 +0100 Subject: [PATCH] Accordion Lattice simple fit --- GaussianBeams_simple_rayleigh.m | 243 ++++++++++++++++++++++++++++++++ 1 file changed, 243 insertions(+) create mode 100644 GaussianBeams_simple_rayleigh.m diff --git a/GaussianBeams_simple_rayleigh.m b/GaussianBeams_simple_rayleigh.m new file mode 100644 index 0000000..7cfa7ea --- /dev/null +++ b/GaussianBeams_simple_rayleigh.m @@ -0,0 +1,243 @@ +%% Elementary constants +clear; +load qf +epsilon=8.854187817e-12; +c = 299792458; %Speed of light [m/s] +h = 6.626070040e-34; %Planck constant [J/Hz] +hbar = h/2/pi; +u0 = 4*pi*1e-7; %Vacuum permeability [N/A^2] +e0 = 1/u0/(c^2); %Vacuum permittivity [A^2*s^4/kg/m^3] +qe = 1.6021766208e-19; %Elementary charge [C] +G = 6.67408e-11; %Gravitational constant [m^3/kg/s^2] +kB = 1.38064852e-23; %Boltzmann constant [J/K] +me = 9.10938356e-31; %Electron rest mass [kg] +mp = 1.672621898e-27; %Proton mass [kg] +uB = qe*hbar /2/me; %Bohr magneton [J/T] +uN = qe*hbar /2/mp; %Nuclear magnton [J/T] +alpha = u0*qe^2*c/(2*h); %Fine structure constant +u = 1.660539040e-27; %Atomic mass unit [kg] +a0 = 4*pi*e0*hbar^2/me/qe^2; %Bohr radius [m] +Eh = hbar^2/(me*a0^2); %Rydberg energy [eV] +%% Laser Properties +wavelength=532e-9; %Laser wavlength [m] +k=2*pi/wavelength; %Wave vector [m^-1] +P=10; %Power of each beam [W] +wz0=60e-6; %Beam waist in z-theta-direction[m] +wy0=250e-6; +zRy=pi*wy0^2./wavelength; %Rayleigh wavelength +zRz=pi*wz0^2./wavelength; +theta=0.5/180*pi; %Half angle of interference +x_ext =wz0./tan(theta); + +%Dysprosium Properties +polar532=350*qe^2*a0^2/Eh; %Dynamical Polarizability of Dy at 532nm +polar1064=193*qe^2*a0^2/Eh; +mDy = 164*u; %Dysprosium mass [kg] +%% Trapping frequency as Plane wave +vz=sqrt(P.*polar532*k^2*sin(theta)^2/(pi^3*wz0*wy0*c*mDy*epsilon)); + +d = wavelength/(2*tan(theta)); %Fringe distance +focus = 100e-3; +D = wavelength*focus/d; +%% Gaussian beams + +%Coordinate System +n=[2001,3,2001]; +xmax=(x_ext.*3/2); +xmin=-(x_ext).*3/2; +ymax=2.*wy0; +ymin=-2.*wy0; +zmax=2.*wz0; +zmin=-2.*wz0; + +for i=1:3 + x=linspace(xmin,xmax,n(1)); + y=linspace(ymin,ymax,n(2)); + z=linspace(zmin,zmax,n(3)); + [X,Y,Z]=meshgrid(x,y,z); + + %Rotation of the coordinte system + xTheta=X.*cos(theta)+Z.*sin(theta); + xPhi=X.*cos(theta)-Z.*sin(theta); + zTheta=-X.*sin(theta)+Z.*cos(theta); + zPhi=X.*sin(theta)+Z.*cos(theta); + + %Electric field + E0=sqrt(2*P/pi/(wz0*wy0)/c/epsilon); %Electric field at 0 + + wz1=wz0*sqrt(1+(xTheta./zRz).^2); %Beam waist in z-direction for Beam 1 + wy1=wy0*sqrt(1+(xTheta./zRy).^2); %Beam waist in y-direction for Beam 1 + wz2=wz0*sqrt(1+(xPhi./zRz).^2); %Beam waist in z-direction for Beam 2 + wy2=wy0*sqrt(1+(xPhi./zRy).^2); %Beam waist in y-direction for Beam 1 + + R1y=xTheta.*(1+(zRy./xTheta).^2); %Rayleigh length Beam 1 + R1z=xTheta.*(1+(zRz./xTheta).^2); + R2y=xPhi.*(1+(zRy./xPhi).^2); %Rayleigh length Beam 2 + R2z=xPhi.*(1+(zRz./xPhi).^2); + + R1y(:,floor(n(1)/2)+1,floor(n(3)/2)+1)=Inf; %Rayleigh length ->0 + R1z(:,floor(n(1)/2)+1,floor(n(3)/2)+1)=Inf; + R2y(:,floor(n(1)/2)+1,floor(n(3)/2)+1)=Inf; + R2z(:,floor(n(1)/2)+1,floor(n(3)/2)+1)=Inf; + + gouyPhase1=atan(wavelength.*xTheta./(pi.*wz0.*wy0)); + gouyPhase2=atan(wavelength.*xPhi./(pi.*wz0.*wy0)); + + + E1=E0.*sqrt(wz0.*wy0)./sqrt(wz1.*wy1).*exp(-(zTheta./wz1).^2-(Y./wy1).^2).*exp(1i*(k.*xTheta-gouyPhase1)+1i.*k.*(zTheta.^2./2./R1z+Y.^2./2./R1y)); %Electric field Beam 1 + E2=E0.*sqrt(wz0.*wy0)./sqrt(wz2.*wy2).*exp(-(zPhi./wz2).^2-(Y./wy2).^2).*exp(1i*(k.*xPhi-gouyPhase2)+1i.*k.*(zPhi.^2./2./R2z+Y.^2./2./R2y)); %Electric field Beam 2 + + + %Intensity and trapping Potential + + Itot=abs(E1+E2).^2./2*c*epsilon; + + if i==1 %Set y=0 + Itot_y=squeeze(Itot(2,:,:)); + Imax=Itot_y(1001,1001); + + figure + imagesc(x,z, Itot_y') + title('x-z Interference Pattern for \theta=', theta*180/pi); + colormap(QF); + caxis([0 Imax]); + + %z-profile for x=0 + Udip=0.5.*Itot_y.*polar532./epsilon./c; + Umid=squeeze(Udip(floor(n(1)/2)+1,:)); + zrange = find((-pi/(3*k*sin(theta))) < z & z < (pi./(3*k*sin(theta)))); + zfit = z(zrange(1):zrange(length(zrange))); + Ufit = Umid(zrange(1):zrange(length(zrange))); + zfit = zfit(:); + Ufit = Ufit(:); + [pks, locs]=findpeaks(Umid); + p = polyfitn(zfit,Ufit,'constant x^2'); + +% figure +% plot(z,Umid, zfit, Ufit,'.', zfit, polyvaln(p,zfit)) +% title('Trapping potential in z-direction for x=y=0') +% %Fitted Trapping frequency +% nu_zfixed=sqrt(-p.Coefficients(2)*2/mDy)/2/pi; + + + %FWHM of x-profile (actually 1/e^2) + Uz0=squeeze(Udip(:,floor(n(1)/2)+1)); + Ufwhm= find(Uz0>Uz0(floor(n(1)/2)+1)/exp(2)); + x_ext_fwhm=abs(x(Ufwhm(1))-x(Ufwhm(length(Ufwhm)))); + x_range=x(Ufwhm); + + %x-profile for z=0 + xrange=find((-x_ext_fwhm/10)