MAJOR CHANGE: Refactored code, added functions to calculate PSD, collision rate

calculateDipoleTrapPotential.py

@@ -66,6 +66,9 @@ def calculateTrapFrequency(w_x, w_z, Power, Polarizability, m = 164*u.u, dir = '
     TrapFrequency = np.nan
     if dir == 'x':
         TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth / (m*w_x**2))).decompose()
+    elif dir == 'y':
+        zReff = np.sqrt(2) * z_R(w_x, 1.064*u.um) * z_R(w_z, 1.064*u.um) / np.sqrt(z_R(w_x, 1.064*u.um)**2 + z_R(w_z, 1.064*u.um)**2)
+        TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(2 * TrapDepth/ (m*zReff**2))).decompose()
     elif dir == 'z':
         TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth/ (m*w_z**2))).decompose()
     return round(TrapFrequency.value, 2)*u.Hz
@@ -84,93 +87,49 @@ def extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, axis):
     dv = pcov[0][0]**0.5 
     return v, dv, popt, pcov
 
-def plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, pcov):
+def meanThermalVelocity(T, m = 164*u.u):
-    v = popt[0]
+    return 4 * np.sqrt((ac.k_B * T) /(np.pi * m))
-    dv = pcov[0][0]**0.5
-    happrox = harmonic_potential(Positions[axis, :].value, *popt)
-    plt.figure()
-    plt.plot(Positions[axis, :].value, happrox, '-r', label = '\u03BD = %.1f \u00B1 %.2f Hz'% tuple([v,dv]))
-    plt.plot(Positions[axis, :], TrappingPotential[axis], 'ob', label = 'Gaussian Potential')
-    plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold')
-    plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold')
-    plt.ylim([-TrapDepthInKelvin.value, max(TrappingPotential[axis].value)])
-    plt.tight_layout()
-    plt.grid(visible=1)
-    plt.legend(prop={'size': 12, 'weight': 'bold'})
-    plt.show()
 
-def plotPotential(Positions, Powers, ComputedPotentials, axis, TrapDepthLabels):
+def particleDensity(w_x, w_z, Power, Polarizability, N, T, m = 164*u.u): # For a thermal cloud
+    v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x')
+    v_y = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'y')
+    v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'z')
 
-    ## plot of the measured parameter vs. scan  parameter
+    return N * (2 * np.pi)**3 * (v_x * v_y * v_z) *  (m / (2 * np.pi * ac.k_B * T))**(3/2)
-    plt.figure(figsize=(9, 7))
-    j = 0
-    for i in range(np.size(ComputedPotentials, 0)):
-        v, dv, popt, pcov = extractTrapFrequency(Positions, ComputedPotentials[i], TrapDepthInKelvin, axis)
-        unit = 'Hz'
-        if v <= 0.0:
-            v        = np.nan
-            dv       = np.nan
-            unit = 'Hz'
-        elif v > 0.0 and orderOfMagnitude(v) > 2:
-            v    = v / 1e3 # in kHz
-            dv   = dv / 1e3 # in kHz
-            unit = 'kHz'
-        tf_label = '\u03BD = %.1f \u00B1 %.2f %s'% tuple([v,dv,unit])
-        if np.size(ComputedPotentials, 0) == len(Powers):
-            plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'P = ' + str(Powers[i]) + ' W; ' + TrapDepthLabels[i] + '; ' + tf_label) 
-        else:
-            if i % 2 == 0 and j < len(Powers):
-                plt.plot(Positions[axis], ComputedPotentials[i][axis], '--',label = 'P = ' + str(Powers[j]) + ' W; ' + TrapDepthLabels[j] + '; ' + tf_label) 
-            elif i % 2 != 0 and j < len(Powers):
-                plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'P = ' + str(Powers[j]) + ' W; ' + tf_label) 
-                j = j + 1
-    if axis == 0:
-        dir = 'X'
-    elif axis == 1:
-        dir = 'Y'
-    else:
-        dir = 'Z'
-    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(loc=3, prop={'size': 12, 'weight': 'bold'})
-    plt.show()
-    # plt.savefig('pot_' + dir + '.png')
 
-if __name__ == '__main__':
+def thermaldeBroglieWavelength(T, m = 164*u.u):
+    return np.sqrt((2*np.pi*ac.hbar**2)/(m*ac.k_B*T))
 
-    # Powers = [0.1, 0.5, 2]
+def scatteringLength(B):
-    # Powers = [5, 20, 40]
+    a_bkg = 87 * ac.a0
-    Powers = [40]
+    #resonanceWidth = 0.005 * u.G
-    Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability
+    #resonanceB = 0.5 * u.G
-    w_x, w_z = 34*u.um, 27.5*u.um # Beam Waists in the x and y directions
-    # w_x, w_z = 30*u.um, 30*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
+    #return a_bkg * (1 - resonanceWidth/(B - resonanceB))
-    extent = 1e4 # range of spatial coordinates in one direction to calculate trap potential over
+    return a_bkg
+
+def dipolarLength(mu = 9.93 * ac.muB, m = 164*u.u):
+    return (m * ac.mu0 * mu**2) / (8 * np.pi * ac.hbar**2)
+
+def scatteringCrossSection(B):
+    return 8 * np.pi * scatteringLength(B)**2 + ((32*np.pi)/45) * dipolarLength()**2
+
+def calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N, T, B): #For a 3D Harmonic Trap
+    return (particleDensity(w_x, w_z, Power, Polarizability, N, T) * scatteringCrossSection(B) * meanThermalVelocity(T) / (2 * np.sqrt(2))).decompose()
+
+def calculatePSD(w_x, w_z, Power, Polarizability, N, T):
+    return (particleDensity(w_x, w_z, Power, Polarizability, N, T, m = 164*u.u) * thermaldeBroglieWavelength(T)**3).decompose()
+
+def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
+
+    axis = options['axis']
+    extent = options['extent']
+    gravity = options['gravity']
+    astigmatism = options['astigmatism']
         
     TrappingPotential = []
-    ComputedPotentials = []
-    TrapDepthLabels = []
-
-    gravity = True
-    astigmatism = True
-
-    tilt_gravity = True
-    theta = 5 # in degrees
-    tilt_axis = [1, 0, 0] # lab space coordinates are rotated about x-axis in reference frame of beam
-
-    disp_foci = 3 * z_R(w_0 = np.asarray([30]), lamb = 1.064)[0]*u.um # difference in position of the foci along the propagation direction (Astigmatism)
-
-    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)
+    IdealTrapDepthInKelvin = (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
 
@@ -189,50 +148,185 @@ if __name__ == '__main__':
     Positions        = np.vstack((x_Positions, y_Positions, z_Positions)) * projection_axis[:, np.newaxis]
 
     IdealTrappingPotential        = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability) 
-        IdealTrappingPotential        = IdealTrappingPotential + np.zeros((3, len(IdealTrappingPotential))) * IdealTrappingPotential.unit
+    IdealTrappingPotential        = IdealTrappingPotential * (np.ones((3, len(IdealTrappingPotential))) * projection_axis[:, np.newaxis])
     IdealTrappingPotential        = (IdealTrappingPotential/ac.k_B).to(u.uK) 
     
     if gravity and not astigmatism:
-            ComputedPotentials.append(IdealTrappingPotential)
         # Influence of Gravity
         m = 164*u.u
         gravity_axis             = np.array([0, 0, 1]) 
+        tilt_gravity             = options['tilt_gravity']
+        theta                    = options['theta']
+        tilt_axis                = options['tilt_axis']
         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        = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability)
+        TrappingPotential        = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
         TrappingPotential        = (TrappingPotential/ac.k_B).to(u.uK)
         
     elif not gravity and astigmatism:
-            ComputedPotentials.append(IdealTrappingPotential)
         # Influence of Astigmatism
+        disp_foci = options['disp_foci']
         TrappingPotential        = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci, alpha = Polarizability)
-            TrappingPotential        = TrappingPotential + np.zeros((3, len(TrappingPotential))) * TrappingPotential.unit
+        TrappingPotential        = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis])
         TrappingPotential        = (TrappingPotential/ac.k_B).to(u.uK)
     
     elif gravity and astigmatism:
-            ComputedPotentials.append(IdealTrappingPotential)
         # Influence of Gravity and Astigmatism
         m = 164*u.u
         gravity_axis             = np.array([0, 0, 1]) 
+        tilt_gravity             = options['tilt_gravity']
+        theta                    = options['theta']
+        tilt_axis                = options['tilt_axis']
+        disp_foci                = options['disp_foci']
         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        = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci, alpha = Polarizability) + gravitational_potential(gravity_axis_positions, m)
+        TrappingPotential        = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci, alpha = Polarizability)
+        TrappingPotential        = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
         TrappingPotential        = (TrappingPotential/ac.k_B).to(u.uK)
     
     else:
         TrappingPotential = IdealTrappingPotential
 
-        # v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x')
+    if TrappingPotential[axis][0] > TrappingPotential[axis][-1]:
-        # v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'z')
+        EffectiveTrapDepthInKelvin = TrappingPotential[axis][-1] - min(TrappingPotential[axis])
+    elif TrappingPotential[axis][0] < TrappingPotential[axis][-1]:
+        EffectiveTrapDepthInKelvin = TrappingPotential[axis][0] - min(TrappingPotential[axis])
 
-        # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, axis)
+    TrapDepthsInKelvin = [IdealTrapDepthInKelvin, EffectiveTrapDepthInKelvin]
-        # plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, pcov)
 
+    v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x')
+    v_y = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'y')
+    v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'z')
+    CalculatedTrapFrequencies = [v_x, v_y, v_z]
+
+    v, dv, popt, pcov = extractTrapFrequency(Positions, IdealTrappingPotential, IdealTrapDepthInKelvin, axis)
+    IdealTrapFrequency = [v, dv]
+    v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, EffectiveTrapDepthInKelvin, axis)
+    TrapFrequency = [v, dv]
+    ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency]
+
+    # Gamma_elastic = calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N = 1.13 * 1e7, T = 22 * u.uK, B = 0 * u.G)
+    # PSD = calculatePSD(w_x, w_z, Power, Polarizability, N = 1.13 * 1e7, T = 22 * u.uK).decompose()
+
+    return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies
+
+def plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, pcov):
+    v = popt[0]
+    dv = pcov[0][0]**0.5
+    happrox = harmonic_potential(Positions[axis, :].value, *popt)
+    plt.figure()
+    plt.plot(Positions[axis, :].value, happrox, '-r', label = '\u03BD = %.1f \u00B1 %.2f Hz'% tuple([v,dv]))
+    plt.plot(Positions[axis, :], TrappingPotential[axis], 'ob', label = 'Gaussian Potential')
+    plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold')
+    plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold')
+    plt.ylim([-TrapDepthInKelvin.value, max(TrappingPotential[axis].value)])
+    plt.tight_layout()
+    plt.grid(visible=1)
+    plt.legend(prop={'size': 12, 'weight': 'bold'})
+    plt.show()
+
+def generate_label(v, dv):
+    unit = 'Hz'
+    if v <= 0.0:
+        v        = np.nan
+        dv       = np.nan
+        unit = 'Hz'
+    elif v > 0.0 and orderOfMagnitude(v) > 2:
+        v    = v / 1e3 # in kHz
+        dv   = dv / 1e3 # in kHz
+        unit = 'kHz'
+    tf_label = '\u03BD = %.1f \u00B1 %.2f %s'% tuple([v,dv,unit])
+    return tf_label
+
+def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterateOver = [], save = False):
+
+    plt.figure(figsize=(9, 7))
+    for i in range(np.size(ComputedPotentials, 0)):
+        
+        if i % 2 == 0:
+            j = int(i / 2)
+        else:
+            j = int((i - 1) / 2)
+        
+        IdealTrapDepthInKelvin = Params[j][0][0]
+        EffectiveTrapDepthInKelvin = Params[j][0][1]
+        
+        idealv = Params[j][2][0][0]
+        idealdv = Params[j][2][0][1]
+        
+        v = Params[j][2][1][0]
+        dv = Params[j][2][1][1]
+        
+        if listToIterateOver:
+            if np.size(ComputedPotentials, 0) == len(listToIterateOver):
+                plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'Trap Depth = ' + str(round(EffectiveTrapDepthInKelvin.value, 2)) + ' ' + str(EffectiveTrapDepthInKelvin.unit) + '; ' + generate_label(v, dv)) 
+            else:
+                if i % 2 == 0:
+                    plt.plot(Positions[axis], ComputedPotentials[i][axis], '--', label = 'Trap Depth = ' + str(round(IdealTrapDepthInKelvin.value, 2)) + ' ' + str(IdealTrapDepthInKelvin.unit) + '; ' + generate_label(idealv, idealdv)) 
+                elif i % 2 != 0:
+                    plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'Effective Trap Depth = ' + str(round(EffectiveTrapDepthInKelvin.value, 2)) + ' ' + str(EffectiveTrapDepthInKelvin.unit) + '; ' + generate_label(v, dv)) 
+        else:
+            if i % 2 == 0:
+                plt.plot(Positions[axis], ComputedPotentials[i][axis], '--', label = 'Trap Depth = ' + str(round(IdealTrapDepthInKelvin.value, 2)) + ' ' + str(IdealTrapDepthInKelvin.unit) + '; ' + generate_label(idealv, idealdv))
+            elif i % 2 != 0:
+                plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'Effective Trap Depth = ' + str(round(EffectiveTrapDepthInKelvin.value, 2)) + ' ' + str(EffectiveTrapDepthInKelvin.unit) + '; ' + generate_label(v, dv))
+    if axis == 0:
+        dir = 'X'
+    elif axis == 1:
+        dir = 'Y'
+    else:
+        dir = 'Z'
+    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(loc=3, prop={'size': 12, 'weight': 'bold'})
+    if save:
+        plt.savefig('pot_' + dir + '.png')
+    plt.show()
+
+if __name__ == '__main__':            
+
+    Power = 40*u.W
+    Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability
+    w_x, w_z = 27.5*u.um, 33.8*u.um # Beam Waists in the x and y directions
+    
+    options = {
+        '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
+        'gravity': True,
+        'astigmatism': False,
+        'tilt_gravity': True,
+        'theta': 5, # in degrees
+        'tilt_axis': [1, 0, 0], # lab space coordinates are rotated about x-axis in reference frame of beam
+        'disp_foci': 3 * z_R(w_0 = np.asarray([30]), lamb = 1.064)[0]*u.um # difference in position of the foci along the propagation direction (Astigmatism)
+    }
+    
+    ComputedPotentials = [] 
+    Params = []  
+
+    Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies = computeTrapPotential(w_x, w_z, Power, Polarizability, options)
+    ComputedPotentials.append(IdealTrappingPotential)
     ComputedPotentials.append(TrappingPotential)
+    Params.append([TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies])
 
     ComputedPotentials = np.asarray(ComputedPotentials)
-    plotPotential(Positions, Powers, ComputedPotentials, axis, TrapDepthLabels)
+    plotPotential(Positions, ComputedPotentials, options['axis'], Params)
+
+    # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, options['axis'])
+    # plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, options['axis'], popt, pcov)
+
+    # Power = [10, 20, 25, 30, 35, 40]*u.W # Single Beam Power
+    # for p in Power: 
+    #     Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies = computeTrapPotential(w_x, w_z, p, Polarizability, options)
+    #     ComputedPotentials.append(IdealTrappingPotential)
+    #     ComputedPotentials.append(TrappingPotential)
+    #     Params.append([TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies])
+        
+    # ComputedPotentials = np.asarray(ComputedPotentials)
+    # plotPotential(Positions, ComputedPotentials, options['axis'], Params)