diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py index 9b3d51e..1ac610e 100644 --- a/calculateDipoleTrapPotential.py +++ b/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): - a_bkg = 87 * ac.a0 - #resonanceWidth = 0.005 * u.G - #resonanceB = 0.5 * u.G +def scatteringLength(B, FR_choice = 1, ABKG_choice = 1): + # Dy 164 a_s versus B in 0 to 8G range + # should match SupMat of PhysRevX.9.021012, fig S5 and descrption + # https://journals.aps.org/prx/supplemental/10.1103/PhysRevX.9.021012/Resubmission_Suppmat.pdf - #return a_bkg * (1 - resonanceWidth/(B - resonanceB)) - return a_bkg + if FR_choice == 1: # new values + + 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('%.2f' % (Gamma_elastic.value) * (Gamma_elastic.unit)) - print('%.2E' % (PSD.value)) + print('Particle Density = %.2E ' % (n.value) + str(n.unit)) + print('Elastic Collision Rate = %.2f ' % (Gamma_elastic.value) + str(Gamma_elastic.unit)) + 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_y) - print(v_z) + print('v_x = %.2f ' %(v_x.value) + str(v_x.unit)) + print('v_y = %.2f ' %(v_y.value) + str(v_y.unit)) + 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)