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.
This commit is contained in:
parent
b9b1321453
commit
4336a657a4
@ -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)
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user