Adding a new script to calculate Dipole Trap characteristics

+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);

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)
+        

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.Power     = 0.2;
-SidebandBeam.Waist     = 16.67e-03;
+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;
+n = 10000;
-t = 2*pi*rand(n,1);
+xmin = -BeamRadius;
-r = BeamRadius*sqrt(rand(n,1));
+xmax = BeamRadius;
-% xmin = -15e-03;
+x    = xmin+rand(n,1)*(xmax-xmin);
-% xmax = 15e-03;
+
-% x    = xmin+rand(n,1)*(xmax-xmin);
+y = ones(n,1) * 0;
-x = r.*cos(t);
+
-y = ones(n,1) * 2e-03;
+zmin = -BeamRadius;
-z = r.*sin(t);
+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;
+rows = 35;
-columns = 100;
+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