From 4336a657a4da84a3131d1d0393b11a3a3a72fa04 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 13 Feb 2023 20:52:01 +0100 Subject: [PATCH] Added correct computation of scattering lengths around Feshbach resonances, added plotting of these scattering lengths vs B field and added the effect due to modulation of beam. --- calculateDipoleTrapPotential.py | 126 ++++++++++++++++++++------------ 1 file changed, 78 insertions(+), 48 deletions(-) diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py index 1ac610e..6f5b020 100644 --- a/calculateDipoleTrapPotential.py +++ b/calculateDipoleTrapPotential.py @@ -36,6 +36,12 @@ def find_nearest(array, value): idx = (np.abs(array - value)).argmin() return idx +def arccos_modulation(mod_amp, n_points): + phi = np.linspace(0, 2*np.pi, n_points) + first_half = 2/np.pi * np.arccos(phi/np.pi-1) - 1 + second_half = np.flip(first_half) + return mod_amp * np.concatenate((first_half, second_half)) + ##################################################################### # BEAM PARAMETERS # ##################################################################### @@ -80,19 +86,10 @@ def scatteringLength(B, FR_choice = 1, ABKG_choice = 1): #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 + 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] * u.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 + 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] * u.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: @@ -104,19 +101,14 @@ def scatteringLength(B, FR_choice = 1, ABKG_choice = 1): #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 + resonanceB = [1.298, 2.802, 3.370, 5.092, 7.154, 2.592, 2.338, 2.177] * u.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 + resonancewB = [0.018, 0.047, 0.048, 0.145, 0.020, 0.008, 0.001, 0.001] * u.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 + #Get scattering length + np.seterr(divide='ignore') + a_s = a_bkg * np.prod([(1 - resonancewB[j] / (B - resonanceB[j])) for j in range(len(resonanceB))]) + return a_s, a_bkg def dipolarLength(mu = 9.93 * ac.muB, m = 164*u.u): return (m * ac.mu0 * mu**2) / (12 * np.pi * ac.hbar**2) @@ -148,15 +140,20 @@ def astigmatic_single_gaussian_beam_potential(positions: "np.ndarray|u.quantity. 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,:] - (del_y/2), waists[0], wavelength))**2 + (positions[2,:]/w(positions[1,:] + (del_y/2), 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 harmonic_potential(pos, v, xoffset, yoffset, m = 164*u.u): U_Harmonic = ((0.5 * m * (2 * np.pi * v*u.Hz)**2 * (pos*u.um - xoffset*u.um)**2)/ac.k_B).to(u.uK) + yoffset*u.uK return U_Harmonic.value +# def modulated_single_gaussian_beam_potential(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", alpha:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: +# mod_amp = 0.5 * waists[0] +# n_points = int(len(positions[0,:])/2) +# x_mod = arccos_modulation(mod_amp, n_points) +# 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.trapz(np.exp(-2 * (np.subtract(x_mod, positions[0,:])/w(positions[1,:], waists[0], wavelength))**2)) +# return U + ##################################################################### # COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS # ##################################################################### @@ -196,7 +193,13 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): extent = options['extent'] gravity = options['gravity'] astigmatism = options['astigmatism'] - + modulation = options['modulation'] + + if modulation: + aspect_ratio = options['aspect_ratio'] + current_ar = w_x/w_z + w_x = w_x * (aspect_ratio / current_ar) + TrappingPotential = [] TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability) IdealTrapDepthInKelvin = (TrapDepth/ac.k_B).to(u.uK) @@ -261,11 +264,13 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): else: TrappingPotential = IdealTrappingPotential - + if TrappingPotential[axis][0] > TrappingPotential[axis][-1]: EffectiveTrapDepthInKelvin = TrappingPotential[axis][-1] - min(TrappingPotential[axis]) elif TrappingPotential[axis][0] < TrappingPotential[axis][-1]: EffectiveTrapDepthInKelvin = TrappingPotential[axis][0] - min(TrappingPotential[axis]) + else: + EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin TrapDepthsInKelvin = [IdealTrapDepthInKelvin, EffectiveTrapDepthInKelvin] @@ -283,7 +288,7 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies ##################################################################### -# PLOT TRAP POTENTIALS # +# PLOTTING # ##################################################################### def plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, axis, popt, pcov): @@ -352,6 +357,8 @@ def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterat dir = 'Y' else: dir = 'Z' + + plt.ylim(top = 0) plt.xlabel(dir + ' Direction (um)', fontsize= 12, fontweight='bold') plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold') plt.tight_layout() @@ -361,8 +368,29 @@ def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterat plt.savefig('pot_' + dir + '.png') plt.show() +def plotScatteringLengths(): + BField = np.arange(0, 2.59, 1e-3) * u.G + a_s_array = np.zeros(len(BField)) * ac.a0 + for idx in range(len(BField)): + a_s_array[idx], a_bkg = scatteringLength(BField[idx]) + rmelmIdx = [i for i, x in enumerate(np.isinf(a_s_array.value)) if x] + for x in rmelmIdx: + a_s_array[x-1] = np.inf * ac.a0 + + plt.figure(figsize=(9, 7)) + plt.plot(BField, a_s_array/ac.a0, '-b') + plt.axhline(y = a_bkg/ac.a0, color = 'r', linestyle = '--') + plt.text(min(BField.value) + 0.5, (a_bkg/ac.a0).value + 1, '$a_{bkg}$ = %.2f a0' %((a_bkg/ac.a0).value), fontsize=14, fontweight='bold') + plt.xlim([min(BField.value), max(BField.value)]) + plt.ylim([65, 125]) + plt.xlabel('B field (G)', fontsize= 12, fontweight='bold') + plt.ylabel('Scattering length (a0)', fontsize= 12, fontweight='bold') + plt.tight_layout() + plt.grid(visible=1) + plt.show() + ##################################################################### -# FUNCTION CALLS BELOW # +# RUN SCRIPT WITH OPTIONS BELOW # ##################################################################### if __name__ == '__main__': @@ -372,12 +400,12 @@ if __name__ == '__main__': 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 + 'axis': 0, # axis referenced to the beam along which you want the dipole trap potential + 'extent': 3e2, # range of spatial coordinates in one direction to calculate trap potential over 'modulation': True, - 'aspect_ratio': 4.6, - 'gravity': True, - 'tilt_gravity': True, + 'aspect_ratio': 3.67, + 'gravity': False, + 'tilt_gravity': False, 'theta': 5, # in degrees 'tilt_axis': [1, 0, 0], # lab space coordinates are rotated about x-axis in reference frame of beam 'astigmatism': False, @@ -387,17 +415,21 @@ if __name__ == '__main__': 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]) + 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, ComputedPotentials, options['axis'], Params) + ComputedPotentials = np.asarray(ComputedPotentials) + plotPotential(Positions, ComputedPotentials, options['axis'], Params) AtomNumber = 1.13 * 1e7 Temperature = 30 * u.uK - BField = 1 * u.G + BField = 2.1 * u.G + + aspect_ratio = 3.67 + init_ar = w_x / w_z + w_x = w_x * (aspect_ratio / init_ar) 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) @@ -415,12 +447,10 @@ if __name__ == '__main__': 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() - + print('a_s = %.2f ' %(scatteringLength(BField)[0] / ac.a0)) + + # plotScatteringLengths() + # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, options['axis']) # plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, options['axis'], popt, pcov)