Cosmetic changes

calculateDipoleTrapPotential.py

@@ -31,6 +31,11 @@ def rotation_matrix(axis, theta):
                      [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                      [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
 
+def find_nearest(array, value):
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return idx
+
 #####################################################################
 # BEAM PARAMETERS                                                   #
 #####################################################################
@@ -44,7 +49,7 @@ def w(pos, w_0, lamb):
     return w_0*np.sqrt(1+(pos / z_R(w_0, lamb))**2)
 
 #####################################################################
-# RELEVANT PARAMETERS FOR EVAPORATIVE COOLING                       #
+# COLLISION RATES, PSD                                              #
 #####################################################################
 
 def meanThermalVelocity(T, m = 164*u.u):
@@ -59,19 +64,65 @@ def particleDensity(w_x, w_z, Power, Polarizability, N, T, m = 164*u.u): # For a
 def thermaldeBroglieWavelength(T, m = 164*u.u):
     return np.sqrt((2*np.pi*ac.hbar**2)/(m*ac.k_B*T))
 
-def scatteringLength(B):
+def scatteringLength(B, FR_choice = 1, ABKG_choice = 1):
-    a_bkg = 87 * ac.a0
+    # Dy 164 a_s versus B in 0 to 8G range
-    #resonanceWidth = 0.005 * u.G
+    # should match SupMat of PhysRevX.9.021012, fig S5 and descrption
-    #resonanceB = 0.5 * u.G
+    # https://journals.aps.org/prx/supplemental/10.1103/PhysRevX.9.021012/Resubmission_Suppmat.pdf
     
-    #return a_bkg * (1 - resonanceWidth/(B - resonanceB))
+    if FR_choice == 1: # new values
-    return a_bkg
+        
+        if ABKG_choice == 1:
+            a_bkg = 85.5 * ac.a0
+        elif ABKG_choice == 2:
+            a_bkg = 93.5 * ac.a0
+        elif ABKG_choice == 3:
+            a_bkg = 77.5 * ac.a0
+        
+        #FR resonances
+        #[B11 B12 B2 B3 B4 B51 B52 B53 B6 B71 B72 B81 B82 B83 B9]
+        resonanceB  = [1.295, 1.306, 2.174, 2.336, 2.591, 2.74, 2.803, 2.78, 3.357, 4.949, 5.083, 7.172, 7.204, 7.134, 76.9] * ac.G #resonance position
+        #[wB11 wB12 wB2 wB3 wB4 wB51 wB52 wB53 wB6 wB71 wB72 wB81 wB82 wB83 wB9]
+        resonancewB = [0.009, 0.010, 0.0005, 0.0005, 0.001, 0.0005, 0.021, 0.015, 0.043, 0.0005, 0.130, 0.024, 0.0005, 0.036, 3.1] * ac.G #resonance width
+        
+        #Get scattering length
+        BField = np.arange(0, 8, 0.5) * ac.G
+        tmp = np.zeros(len(resonanceB)) * ac.a0
+        for idx in range(len(resonanceB)):
+            tmp[idx] = [(1 - resonancewB[idx] / (BField[j] - resonanceB[idx])) for j in range(len(BField))]
+        a_s_array = tmp
+        #index = find_nearest(BField.value, B.value)
+        a_s = 1 #a_s_array[index]
+
+    else: # old values
+        
+        if ABKG_choice == 1:
+            a_bkg = 87.2 * ac.a0
+        elif ABKG_choice == 2:
+            a_bkg = 95.2 * ac.a0
+        elif ABKG_choice == 3:
+            a_bkg = 79.2 * ac.a0
+
+        #FR resonances
+        #[B1 B2 B3 B4 B5 B6 B7 B8]
+        resonanceB  = [1.298, 2.802, 3.370, 5.092, 7.154, 2.592, 2.338, 2.177] * ac.G #resonance position
+        #[wB1 wB2 wB3 wB4 wB5 wB6 wB7 wB8]
+        resonancewB = [0.018, 0.047, 0.048, 0.145, 0.020, 0.008, 0.001, 0.001] * ac.G #resonance width
+        
+        #Get scattering length
+        BField = np.arange(0,8, 0.0001) * ac.G
+        a_s_array = np.zeros(len(BField)) * ac.a0
+        for idx in range(len(BField)):
+            a_s_array[idx] = a_bkg * (1 - resonancewB[idx] / (BField[idx] - resonanceB[idx]))
+        index = find_nearest(BField.value, B.value)
+        a_s = a_s_array[index]
+
+    return a_s, a_s_array, BField
 
 def dipolarLength(mu = 9.93 * ac.muB, m = 164*u.u):
     return (m * ac.mu0 * mu**2) / (12 * np.pi * ac.hbar**2)
 
 def scatteringCrossSection(B):
-    return 8 * np.pi * scatteringLength(B)**2 + ((32*np.pi)/45) * dipolarLength()**2
+    return 8 * np.pi * scatteringLength(B)[0]**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()
@@ -320,17 +371,16 @@ if __name__ == '__main__':
     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
     
-    AspectRatio = 4.6
-    w_x = AspectRatio * w_x
-
     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
+        'modulation': True,
+        'aspect_ratio': 4.6,
         '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
+        'astigmatism': False,
         '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)
     }
     
@@ -347,23 +397,29 @@ if __name__ == '__main__':
 
     AtomNumber = 1.13 * 1e7
     Temperature = 30 * u.uK
-    BField = 0 * u.G
+    BField = 1 * u.G
 
     n = particleDensity(w_x, w_z, Power, Polarizability, N = AtomNumber, T = Temperature, m = 164*u.u).decompose().to(u.cm**(-3))
     Gamma_elastic = calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N = AtomNumber, T = Temperature, B = BField)
     PSD = calculatePSD(w_x, w_z, Power, Polarizability, N = AtomNumber, T = Temperature).decompose()
     
-    print('%.2E' % (n.value)* (n.unit))
+    print('Particle Density = %.2E ' % (n.value) + str(n.unit))
-    print('%.2f' % (Gamma_elastic.value) * (Gamma_elastic.unit))
+    print('Elastic Collision Rate = %.2f ' % (Gamma_elastic.value) + str(Gamma_elastic.unit))
-    print('%.2E' % (PSD.value))
+    print('PSD = %.2E ' % (PSD.value))
 
     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')
 
-    print(v_x)
+    print('v_x = %.2f ' %(v_x.value) + str(v_x.unit))
-    print(v_y)
+    print('v_y = %.2f ' %(v_y.value) + str(v_y.unit))
-    print(v_z)
+    print('v_z = %.2f ' %(v_z.value) + str(v_z.unit))
+
+    #plt.figure()
+    ret = scatteringLength(1 * ac.G)
+    print(ret[1])
+    #plt.plot(ret[2], ret[1])
+    #plt.show()
 
     # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, options['axis'])
     # plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, options['axis'], popt, pcov)