From 435fbccbfa04326f63e00e961e621d58fbfb3114 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 11 Jan 2023 18:54:23 +0100 Subject: [PATCH] Adding a new script to calculate Dipole Trap characteristics --- +Scripts/optimizingForSidebandEnhancement.m | 6 +- calculateDipoleTrapPotential.py | 151 ++++++++++++++++++++ test_MOTCaptureProcessSimulation.m | 43 +++--- 3 files changed, 180 insertions(+), 20 deletions(-) create mode 100644 calculateDipoleTrapPotential.py diff --git a/+Scripts/optimizingForSidebandEnhancement.m b/+Scripts/optimizingForSidebandEnhancement.m index 0639587..4af4332 100644 --- a/+Scripts/optimizingForSidebandEnhancement.m +++ b/+Scripts/optimizingForSidebandEnhancement.m @@ -173,7 +173,7 @@ OptionsStruct.RescalingFactorForSecondParameter = 1000; OptionsStruct.YLabelString = 'Beam Waist (mW)'; OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-10; % OptionsStruct.ZLabelString = 'Enhancement Factor (\eta)'; -OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)'; +OptionsStruct.ZLabelString = 'Loading rate (x 10^{10} atoms/s)'; OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100); options = Helper.convertstruct2cell(OptionsStruct); @@ -189,7 +189,7 @@ SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideban SidebandBeam.Power = 0.2; NumberOfPointsForFirstParam = 20; %iterations of the simulation NumberOfPointsForSecondParam = 20; -DetuningArray = linspace(-1.0, -5.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; +DetuningArray = linspace(-0.5, -2.5, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; BeamWaistArray = linspace(10, 25, NumberOfPointsForSecondParam) * 1e-03; tStart = tic; @@ -229,7 +229,7 @@ OptionsStruct.RescalingFactorForSecondParameter = 1000; OptionsStruct.YLabelString = 'Beam Waist (mW)'; OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-10; % OptionsStruct.ZLabelString = 'Enhancement Factor (\eta)'; -OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)'; +OptionsStruct.ZLabelString = 'Loading rate (x 10^{10} atoms/s)'; OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100); options = Helper.convertstruct2cell(OptionsStruct); diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py new file mode 100644 index 0000000..1e26835 --- /dev/null +++ b/calculateDipoleTrapPotential.py @@ -0,0 +1,151 @@ +import math +import numpy as np +import matplotlib.pyplot as plt +from astropy import units as u, constants as ac + +def rotation_matrix(axis, theta): + """ + Return the rotation matrix associated with counterclockwise rotation about + the given axis by theta radians. + + In 2-D it is just, + thetaInRadians = np.radians(theta) + c, s = np.cos(thetaInRadians), np.sin(thetaInRadians) + R = np.array(((c, -s), (s, c))) + + In 3-D, one way to do it is use the Euler-Rodrigues Formula as is done here + """ + axis = np.asarray(axis) + axis = axis / math.sqrt(np.dot(axis, axis)) + a = math.cos(theta / 2.0) + b, c, d = -axis * math.sin(theta / 2.0) + aa, bb, cc, dd = a * a, b * b, c * c, d * d + bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d + return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)], + [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)], + [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) + +# Rayleigh range +def z_R(w_0:np.ndarray, lamb:float)->np.ndarray: + return np.pi*w_0**2/lamb + +# Beam Radius +def w(pos, w_0, lamb): + return w_0*np.sqrt(1+(pos*lamb/(np.pi*w_0**2))**2) + +def trap_depth(w_1:"float|u.quantity.Quantity", w_2:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity", alpha:float)->"float|u.quantity.Quantity": + return 2*P/(np.pi*w_1*w_2) * (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) + +def gravitational_potential(positions: "np.ndarray|u.quantity.Quantity", m:"float|u.quantity.Quantity"): + return m * ac.g0 * positions + +def single_gaussian_beam_potential(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", P:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um, alpha:"float|u.quantity.Quantity"=184.4)->np.ndarray: + A = 2*P/(np.pi*w(positions[1,:], waists[0], wavelength)*w(positions[1,:], waists[1], wavelength)) + U_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) + U = - U_tilde * A * np.exp(-2 * ((positions[0,:]/w(positions[1,:], waists[0], wavelength))**2 + (positions[2,:]/w(positions[1,:], waists[1], wavelength))**2)) + return U + +def single_gaussian_beam_potential_harmonic_approximation(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", depth:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: + U = - depth * (1 - (2 * (positions[0,:]/waists[0])**2) - (2 * (positions[2,:]/waists[1])**2) - (0.5 * positions[1,:]**2 * np.sum(np.reciprocal(z_R(waists, wavelength)))**2)) + return U + +def plotPotential(Positions, Powers, ComputedPotentials, axis, TrapDepthLabels): + + ## plot of the measured parameter vs. scan parameter + plt.figure(figsize=(9, 7)) + for i in range(np.size(ComputedPotentials, 0)): + plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'P = ' + str(Powers[i]) + ' W; ' + TrapDepthLabels[i]) + if axis == 0: + dir = 'X' + elif axis == 1: + dir = 'Y' + else: + dir = 'Z' + # maxPotentialValue = max(ComputedPotentials.flatten()) + # minPotentialValue = min(ComputedPotentials.flatten()) + # PotentialValueRange = maxPotentialValue - minPotentialValue + # upperlimit = 5 + # if maxPotentialValue > 0: + # upperlimit = maxPotentialValue + # lowerlimit = min(ComputedPotentials.flatten()) - PotentialValueRange/6 + # plt.ylim([lowerlimit, upperlimit]) + plt.xlabel(dir + ' Direction (um)', fontsize= 12, fontweight='bold') + plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold') + plt.tight_layout() + plt.grid(visible=1) + plt.legend(prop={'size': 12, 'weight': 'bold'}) + plt.tight_layout() + # plt.show() + plt.savefig('pot_' + dir + '.png') + +if __name__ == '__main__': + + # Powers = [0.1, 0.5, 2] + # Powers = [5, 10, 20, 30, 40] + Powers = [40] + Polarizability = 160 # in a.u., should we use alpha = 136 or 160 or 184.4? + # w_x, w_z = 34*u.um, 27.5*u.um # Beam Waists in the x and y directions + w_x, w_z = 35*u.um, 35*u.um # Beam Waists in the x and y directions + # w_x, w_z = 20.5*u.um, 20.5*u.um + + axis = 1 # axis referenced to the beam along which you want the dipole trap potential + extent = 1e4 # range of spatial coordinates in one direction to calculate trap potential over + + TrappingPotential = [] + ComputedPotentials = [] + TrapDepthLabels = [] + + gravity = False + astigmatism = False + + tilt_gravity = True + theta = 0 # in degrees + tilt_axis = [1, 0, 0] # lab space coordinates are rotated about x-axis in reference frame of beam + + for p in Powers: + + Power = p*u.W # Single Beam Power + + TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability) + TrapDepthInKelvin = (TrapDepth/ac.k_B).to(u.uK) + TrapDepthLabels.append("Trap Depth = " + str(round(TrapDepthInKelvin.value, 2)) + " " + str(TrapDepthInKelvin.unit)) + + projection_axis = np.array([0, 1, 0]) # default + if axis == 0: + projection_axis = np.array([1, 0, 0]) # radial direction (X-axis) + elif axis == 1: + projection_axis = np.array([0, 1, 0]) # propagation direction (Y-axis) + elif axis == 2: + projection_axis = np.array([0, 0, 1]) # vertical direction (Z-axis) + + x_Positions = np.arange(-extent, extent, 1)*u.um + y_Positions = np.arange(-extent, extent, 1)*u.um + z_Positions = np.arange(-extent, extent, 1)*u.um + Positions = np.vstack((x_Positions, y_Positions, z_Positions)) * projection_axis[:, np.newaxis] + + if not gravity and not astigmatism: + TrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability) + TrappingPotential = TrappingPotential + np.zeros((3, len(TrappingPotential))) * TrappingPotential.unit + TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) + + elif gravity and not astigmatism: + # Influence of Gravity + m = 164*u.u + gravity_axis = np.array([0, 0, 1]) + if tilt_gravity: + R = rotation_matrix(tilt_axis, np.radians(theta)) + gravity_axis = np.dot(R, gravity_axis) + gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis] + TrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability) + gravitational_potential(gravity_axis_positions, m) + TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) + + ComputedPotentials.append(TrappingPotential) + + ComputedPotentials = np.asarray(ComputedPotentials) + plotPotential(Positions, Powers, ComputedPotentials, axis, TrapDepthLabels) + +# Influence of Astigmatism + +# TrappingPotential = single_gaussian_beam_potential_harmonic_approximation(Positions, np.asarray([w_x.value, w_z.value])*u.um, depth = TrapDepth) +# TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) + \ No newline at end of file diff --git a/test_MOTCaptureProcessSimulation.m b/test_MOTCaptureProcessSimulation.m index e77b818..2c027ad 100644 --- a/test_MOTCaptureProcessSimulation.m +++ b/test_MOTCaptureProcessSimulation.m @@ -25,16 +25,16 @@ MOT2D = Simulator.TwoDimensionalMOT(options{:}); Beams = MOT2D.Beams; %% - Run Simulation -MOT2D.NumberOfAtoms = 5000; +MOT2D.NumberOfAtoms = 10000; MOT2D.SidebandBeam = true; MOT2D.PushBeam = false; CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; CoolingBeam.Power = 0.2; -CoolingBeam.Waist = 16.67e-03; +CoolingBeam.Waist = 20e-03; CoolingBeam.Detuning = -1.33*Helper.PhysicsConstants.BlueLinewidth; SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)}; -SidebandBeam.Power = 0.4; -SidebandBeam.Waist = 16.67e-03; +SidebandBeam.Power = 0.2; +SidebandBeam.Waist = 20e-03; SidebandBeam.Detuning = -2.66*Helper.PhysicsConstants.BlueLinewidth; PushBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), Beams)}; PushBeam.Power = 0.025; @@ -227,7 +227,7 @@ clear OptionsStruct %% COOLING BEAM WAIST VS DETUNING MOT2D.NumberOfAtoms = 20000; -MOT2D.MagneticGradient = 0.38; +MOT2D.MagneticGradient = 0.40; MOT2D.SidebandBeam = false; MOT2D.PushBeam = false; CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; @@ -320,15 +320,24 @@ CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; SaturationIntensity = CoolingBeam.SaturationIntensity; SaturationParameter = 0.1 * (8 * Power) / (pi*BeamWaist^2 * SaturationIntensity); % two beams are reflected -n = 1000; -t = 2*pi*rand(n,1); -r = BeamRadius*sqrt(rand(n,1)); -% xmin = -15e-03; -% xmax = 15e-03; -% x = xmin+rand(n,1)*(xmax-xmin); -x = r.*cos(t); -y = ones(n,1) * 2e-03; -z = r.*sin(t); +n = 10000; +xmin = -BeamRadius; +xmax = BeamRadius; +x = xmin+rand(n,1)*(xmax-xmin); + +y = ones(n,1) * 0; + +zmin = -BeamRadius; +zmax = BeamRadius; +z = zmin+rand(n,1)*(zmax-zmin); + + +% t = 2*pi*rand(n,1); +% r = BeamRadius*sqrt(rand(n,1)); +% x = r.*cos(t); +% y = ones(n,1) * 0; +% z = r.*sin(t); + PositionVector = horzcat(x, y, z); %scatter3(zeros(n,1), y, z) CoolingBeamLocalSaturationIntensity = @(x) MOT2D.calculateLocalSaturationIntensity(0.25 * SaturationParameter, x, Origin, WaveVectorEndPoint(BeamNumber,:), BeamRadius, BeamWaist); IntensityProfile = zeros(n,1); @@ -337,8 +346,8 @@ for i=1:n end v = IntensityProfile; % Extract intensity value -rows = 100; -columns = 100; +rows = 35; +columns = 35; Image = zeros(rows, columns); for k = 1 : length(x) row = ceil((x(k) - min(x)) * columns / (max(x) - min(x))); @@ -395,7 +404,7 @@ f_h.Units = 'pixels'; set(0,'units','pixels'); screensize = get(0,'ScreenSize'); f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600]; -[xq,zq] = meshgrid(linspace(-BeamWaist, BeamWaist, n), linspace(-BeamWaist, BeamWaist, n)); +[xq,zq] = meshgrid(linspace(-BeamRadius, BeamRadius, n), linspace(-BeamRadius, BeamRadius, n)); vq = griddata(x,z,v,xq,zq); mesh(xq,zq,vq) hold on