karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Corrected minor bugs, added detailed plotting of crossed beam potential in different directions.

Browse Source
This commit is contained in:
Karthik 2023-03-15 11:29:21 +01:00
parent 6f018c0281
commit 6976e0ef91
2 changed files with 443 additions and 132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

417
ODTCalculations.ipynb
View File

File diff suppressed because one or more lines are too long

158
calculateDipoleTrapPotential.py
View File

@ -6,6 +6,10 @@ from scipy import signal, interpolate
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
from astropy import units as u, constants as ac from astropy import units as u, constants as ac
DY_POLARIZABILITY = 184.4 # in a.u, most precise measured value of Dy polarizability
DY_MASS = 164*u.u
DY_DIPOLE_MOMENT = 9.93 * ac.muB
##################################################################### #####################################################################
# HELPER FUNCTIONS # # HELPER FUNCTIONS #
##################################################################### #####################################################################
@ -80,20 +84,20 @@ def w(pos, w_0, lamb):
# COLLISION RATES, PSD # # COLLISION RATES, PSD #
##################################################################### #####################################################################
def calculateAtomNumber(NCount, pixel_size = 5.86 * u.um, magnification = 0.5, eta = 0.5): def calculateAtomNumber(NCount, pixel_size, magnification, eta):
sigma = 8.468e-14 * (u.m)**(2) sigma = 8.468e-14 * (u.m)**(2)
return (1/eta * 1/sigma * NCount * pixel_size**2/magnification**2).decompose() return (1/eta * 1/sigma * NCount * pixel_size**2/magnification**2).decompose()
def meanThermalVelocity(T, m = 164*u.u): def meanThermalVelocity(T, m = DY_MASS):
return 4 * np.sqrt((ac.k_B * T) /(np.pi * m)) return 4 * np.sqrt((ac.k_B * T) /(np.pi * m))
def particleDensity(w_x, w_z, Power, N, T, m = 164*u.u): # For a thermal cloud def particleDensity(w_x, w_z, Power, N, T, m = DY_MASS): # For a thermal cloud
v_x = calculateTrapFrequency(w_x, w_z, Power, dir = 'x') v_x = calculateTrapFrequency(w_x, w_z, Power, dir = 'x')
v_y = calculateTrapFrequency(w_x, w_z, Power, dir = 'y') v_y = calculateTrapFrequency(w_x, w_z, Power, dir = 'y')
v_z = calculateTrapFrequency(w_x, w_z, Power, dir = 'z') v_z = calculateTrapFrequency(w_x, w_z, Power, dir = 'z')
return N * (2 * np.pi)**3 * (v_x * v_y * v_z) * (m / (2 * np.pi * ac.k_B * T))**(3/2) return N * (2 * np.pi)**3 * (v_x * v_y * v_z) * (m / (2 * np.pi * ac.k_B * T))**(3/2)
def calculateParticleDensityFromMeasurements(v_x, dv_x, v_y, dv_y, v_z, dv_z, w_x, w_z, T_x, T_y, dT_x, dT_y, modulation_depth, N, m = 164*u.u): def calculateParticleDensityFromMeasurements(v_x, dv_x, v_y, dv_y, v_z, dv_z, w_x, w_z, T_x, T_y, dT_x, dT_y, modulation_depth, N, m = DY_MASS):
alpha_x = [(v_x[0]/x)**(2/3) for x in v_x] alpha_x = [(v_x[0]/x)**(2/3) for x in v_x]
dalpha_x = [alpha_x[i] * np.sqrt((dv_x[0]/v_x[0])**2 + (dv_x[i]/v_x[i])**2) for i in range(len(v_x))] dalpha_x = [alpha_x[i] * np.sqrt((dv_x[0]/v_x[0])**2 + (dv_x[i]/v_x[i])**2) for i in range(len(v_x))]
alpha_y = [(v_z[0]/y)**2 for y in v_z] alpha_y = [(v_z[0]/y)**2 for y in v_z]
@ -118,7 +122,7 @@ def calculateParticleDensityFromMeasurements(v_x, dv_x, v_y, dv_y, v_z, dv_z, w_
return pd, dpd, avg_T, avg_dT, new_aspect_ratio return pd, dpd, avg_T, avg_dT, new_aspect_ratio
def thermaldeBroglieWavelength(T, m = 164*u.u): def thermaldeBroglieWavelength(T, m = DY_MASS):
return np.sqrt((2*np.pi*ac.hbar**2)/(m*ac.k_B*T)) return np.sqrt((2*np.pi*ac.hbar**2)/(m*ac.k_B*T))
def scatteringLength(B, FR_choice = 1, ABKG_choice = 1): def scatteringLength(B, FR_choice = 1, ABKG_choice = 1):
@ -161,7 +165,7 @@ def scatteringLength(B, FR_choice = 1, ABKG_choice = 1):
a_s = a_bkg * np.prod([(1 - resonancewB[j] / (B - resonanceB[j])) for j in range(len(resonanceB))]) a_s = a_bkg * np.prod([(1 - resonancewB[j] / (B - resonanceB[j])) for j in range(len(resonanceB))])
return a_s, a_bkg return a_s, a_bkg
def dipolarLength(mu = 9.93 * ac.muB, m = 164*u.u): def dipolarLength(mu = DY_DIPOLE_MOMENT, m = DY_MASS):
return (m * ac.mu0 * mu**2) / (12 * np.pi * ac.hbar**2) return (m * ac.mu0 * mu**2) / (12 * np.pi * ac.hbar**2)
def scatteringCrossSection(B): def scatteringCrossSection(B):
@ -213,19 +217,19 @@ def convert_modulation_depth_to_temperature(modulation_depth):
def gravitational_potential(positions, m): def gravitational_potential(positions, m):
return m * ac.g0 * positions return m * ac.g0 * positions
def single_gaussian_beam_potential(positions, waists, alpha = 184.4, P=1, wavelength=1.064*u.um): def single_gaussian_beam_potential(positions, waists, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um):
A = 2*P/(np.pi*w(positions[1,:], waists[0], wavelength)*w(positions[1,:], waists[1], wavelength)) 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_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,:], waists[0], wavelength))**2 + (positions[2,:]/w(positions[1,:], waists[1], wavelength))**2)) U = - U_tilde * A * np.exp(-2 * ((positions[0,:]/w(positions[1,:], waists[0], wavelength))**2 + (positions[2,:]/w(positions[1,:], waists[1], wavelength))**2))
return U return U
def astigmatic_single_gaussian_beam_potential(positions, waists, del_y, alpha = 184.4, P=1, wavelength=1.064*u.um): def astigmatic_single_gaussian_beam_potential(positions, waists, del_y, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um):
A = 2*P/(np.pi*w(positions[1,:] - (del_y/2), waists[0], wavelength)*w(positions[1,:] + (del_y/2), waists[1], wavelength)) A = 2*P/(np.pi*w(positions[1,:] - (del_y/2), waists[0], wavelength)*w(positions[1,:] + (del_y/2), waists[1], wavelength))
U_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) 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)) 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 return U
def modulated_single_gaussian_beam_potential(positions, waists, alpha = 184.4, P=1, wavelength=1.064*u.um, mod_amp=1): def modulated_single_gaussian_beam_potential(positions, waists, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um, mod_amp=1):
mod_amp = mod_amp * waists[0] mod_amp = mod_amp * waists[0]
n_points = len(positions[0,:]) n_points = len(positions[0,:])
dx, x_mod = modulation_function(mod_amp, n_points, func = 'arccos') dx, x_mod = modulation_function(mod_amp, n_points, func = 'arccos')
@ -237,7 +241,7 @@ def modulated_single_gaussian_beam_potential(positions, waists, alpha = 184.4, P
U = - U_tilde * A * 1/(2*mod_amp) * np.trapz(dU, dx = dx, axis = 0) U = - U_tilde * A * 1/(2*mod_amp) * np.trapz(dU, dx = dx, axis = 0)
return U return U
def harmonic_potential(pos, v, xoffset, yoffset, m = 164*u.u): def harmonic_potential(pos, v, xoffset, yoffset, m = DY_MASS):
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 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 return U_Harmonic.value
@ -245,38 +249,60 @@ def gaussian_potential(pos, amp, waist, xoffset, yoffset):
U_Gaussian = amp * np.exp(-2 * ((pos + xoffset) / waist)**2) + yoffset U_Gaussian = amp * np.exp(-2 * ((pos + xoffset) / waist)**2) + yoffset
return U_Gaussian return U_Gaussian
def crossed_beam_potential(positions, theta, waists, P, alpha = 184.4, wavelength=1.064*u.um): def crossed_beam_potential(positions, waists, P, options, alpha = DY_POLARIZABILITY, wavelength=1.064*u.um):
beam_1_positions = positions delta = options['delta']
A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:], waists[0][0], wavelength)*w(beam_1_positions[1,:], waists[0][1], wavelength))
foci_shift = options['foci_shift']
focus_shift_beam_1 = foci_shift[0]
focus_shift_beam_2 = foci_shift[1]
beam_disp = options['beam_disp']
beam_1_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[0])).T * beam_disp[0].unit
beam_2_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[1])).T * beam_disp[1].unit
beam_1_positions = positions + beam_1_disp
A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][0], wavelength)*w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][1], wavelength))
U_1_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) U_1_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3)
U_1 = - U_1_tilde * A_1 * np.exp(-2 * ((beam_1_positions[0,:]/w(beam_1_positions[1,:], waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:], waists[0][1], wavelength))**2)) U_1 = - U_1_tilde * A_1 * np.exp(-2 * ((beam_1_positions[0,:]/w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][1], wavelength))**2))
R = rotation_matrix([0, 0, 1], np.radians(theta)) R = rotation_matrix([0, 0, 1], np.radians(delta))
beam_2_positions = np.dot(R, beam_1_positions) beam_2_positions = np.dot(R, beam_1_positions) + beam_2_disp
A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:], waists[1][0], wavelength)*w(beam_2_positions[1,:], waists[1][1], wavelength)) A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][0], wavelength)*w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][1], wavelength))
U_2_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) U_2_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3)
U_2 = - U_2_tilde * A_2 * np.exp(-2 * ((beam_2_positions[0,:]/w(beam_2_positions[1,:], waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:], waists[1][1], wavelength))**2)) U_2 = - U_2_tilde * A_2 * np.exp(-2 * ((beam_2_positions[0,:]/w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][1], wavelength))**2))
U = U_1 + U_2 U = U_1 + U_2
return U return U
def astigmatic_crossed_beam_potential(positions, theta, waists, P, del_y, alpha = 184.4, wavelength=1.064*u.um): def astigmatic_crossed_beam_potential(positions, waists, P, options, alpha = DY_POLARIZABILITY, wavelength=1.064*u.um):
delta = options['delta']
del_y = options['foci_disp_crossed']
del_y_1 = del_y[0] del_y_1 = del_y[0]
del_y_2 = del_y[1] del_y_2 = del_y[1]
beam_1_positions = positions
A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] - (del_y_1/2), waists[0][0], wavelength)*w(beam_1_positions[1,:] + (del_y_1/2), waists[0][1], wavelength))
U_1_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3)
U_1 = - U_1_tilde * A_1 * np.exp(-2 * ((beam_1_positions[0,:]/w(beam_1_positions[1,:] - (del_y_1/2), waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + (del_y_1/2), waists[0][1], wavelength))**2))
R = rotation_matrix([0, 0, 1], np.radians(theta)) foci_shift = options['foci_shift']
beam_2_positions = np.dot(R, beam_1_positions) focus_shift_beam_1 = foci_shift[0]
A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:] - (del_y_2/2), waists[1][0], wavelength)*w(beam_2_positions[1,:] + (del_y_2/2), waists[1][1], wavelength)) focus_shift_beam_2 = foci_shift[1]
beam_disp = options['beam_disp']
beam_1_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[0])).T * beam_disp[0].unit
beam_2_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[1])).T * beam_disp[1].unit
beam_1_positions = positions + beam_1_disp
A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] - (del_y_1/2) + focus_shift_beam_1, waists[0][0], wavelength)*w(beam_1_positions[1,:] + (del_y_1/2) + focus_shift_beam_1, waists[0][1], wavelength))
U_1_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3)
U_1 = - U_1_tilde * A_1 * np.exp(-2 * ((beam_1_positions[0,:]/w(beam_1_positions[1,:] - (del_y_1/2) + focus_shift_beam_1, waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + (del_y_1/2) + focus_shift_beam_1, waists[0][1], wavelength))**2))
R = rotation_matrix([0, 0, 1], np.radians(delta))
beam_2_positions = np.dot(R, beam_1_positions) + beam_2_disp
A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:] - (del_y_2/2) + focus_shift_beam_2, waists[1][0], wavelength)*w(beam_2_positions[1,:] + (del_y_2/2) + focus_shift_beam_2, waists[1][1], wavelength))
U_2_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) U_2_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3)
U_2 = - U_2_tilde * A_2 * np.exp(-2 * ((beam_2_positions[0,:]/w(beam_2_positions[1,:] - (del_y_2/2), waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + (del_y_2/2), waists[1][1], wavelength))**2)) U_2 = - U_2_tilde * A_2 * np.exp(-2 * ((beam_2_positions[0,:]/w(beam_2_positions[1,:] - (del_y_2/2) + focus_shift_beam_2, waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + (del_y_2/2) + focus_shift_beam_2, waists[1][1], wavelength))**2))
U = U_1 + U_2 U = U_1 + U_2
@ -286,10 +312,10 @@ def astigmatic_crossed_beam_potential(positions, theta, waists, P, del_y, alpha
# COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS # # COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS #
##################################################################### #####################################################################
def trap_depth(w_1, w_2, P, alpha = 184.4): def trap_depth(w_1, w_2, P, alpha = DY_POLARIZABILITY):
return 2*P/(np.pi*w_1*w_2) * (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) return 2*P/(np.pi*w_1*w_2) * (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3)
def calculateTrapFrequency(w_x, w_z, Power, dir = 'x', m = 164*u.u): def calculateTrapFrequency(w_x, w_z, Power, dir = 'x', m = DY_MASS):
TrapDepth = trap_depth(w_x, w_z, Power) TrapDepth = trap_depth(w_x, w_z, Power)
TrapFrequency = np.nan TrapFrequency = np.nan
if dir == 'x': if dir == 'x':
@ -301,7 +327,7 @@ def calculateTrapFrequency(w_x, w_z, Power, dir = 'x', m = 164*u.u):
TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth/ (m*w_z**2))).decompose() TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth/ (m*w_z**2))).decompose()
return round(TrapFrequency.value, 2)*u.Hz return round(TrapFrequency.value, 2)*u.Hz
def calculateCrossedBeamTrapFrequency(delta, Waists, Powers, dir = 'x', m = 164*u.u, wavelength=1.064*u.um): def calculateCrossedBeamTrapFrequency(delta, Waists, Powers, dir = 'x', m = DY_MASS, wavelength=1.064*u.um):
TrapDepth_1 = trap_depth(Waists[0][0], Waists[1][0], Powers[0]) TrapDepth_1 = trap_depth(Waists[0][0], Waists[1][0], Powers[0])
TrapDepth_2 = trap_depth(Waists[0][1], Waists[1][1], Powers[1]) TrapDepth_2 = trap_depth(Waists[0][1], Waists[1][1], Powers[1])
@ -392,7 +418,7 @@ def computeTrapPotential(w_x, w_z, Power, options):
if gravity and not astigmatism: if gravity and not astigmatism:
# Influence of Gravity # Influence of Gravity
m = 164*u.u m = DY_MASS
gravity_axis = np.array([0, 0, 1]) gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity'] tilt_gravity = options['tilt_gravity']
theta = options['theta'] theta = options['theta']
@ -407,24 +433,24 @@ def computeTrapPotential(w_x, w_z, Power, options):
elif not gravity and astigmatism: elif not gravity and astigmatism:
# Influence of Astigmatism # Influence of Astigmatism
disp_foci = options['disp_foci'] foci_disp_single = options['foci_disp_single']
TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci) TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = foci_disp_single)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis])
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif gravity and astigmatism: elif gravity and astigmatism:
# Influence of Gravity and Astigmatism # Influence of Gravity and Astigmatism
m = 164*u.u m = DY_MASS
gravity_axis = np.array([0, 0, 1]) gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity'] tilt_gravity = options['tilt_gravity']
theta = options['theta'] theta = options['theta']
tilt_axis = options['tilt_axis'] tilt_axis = options['tilt_axis']
disp_foci = options['disp_foci'] foci_disp_single = options['foci_disp_single']
if tilt_gravity: if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta)) R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis) gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis] 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) TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = foci_disp_single)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m) 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) TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
@ -472,15 +498,14 @@ def computeTrapPotential(w_x, w_z, Power, options):
return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies
else: else:
delta = options['delta']
waists = np.vstack((np.asarray([w_x[0].value, w_z[0].value])*u.um, np.asarray([w_x[1].value, w_z[1].value])*u.um)) waists = np.vstack((np.asarray([w_x[0].value, w_z[0].value])*u.um, np.asarray([w_x[1].value, w_z[1].value])*u.um))
IdealTrappingPotential = crossed_beam_potential(Positions, delta, waists, P = Power) IdealTrappingPotential = crossed_beam_potential(Positions, waists, Power, options)
IdealTrappingPotential = IdealTrappingPotential * (np.ones((3, len(IdealTrappingPotential))) * projection_axis[:, np.newaxis]) IdealTrappingPotential = IdealTrappingPotential * (np.ones((3, len(IdealTrappingPotential))) * projection_axis[:, np.newaxis])
IdealTrappingPotential = (IdealTrappingPotential/ac.k_B).to(u.uK) IdealTrappingPotential = (IdealTrappingPotential/ac.k_B).to(u.uK)
if gravity and not astigmatism: if gravity and not astigmatism:
# Influence of Gravity # Influence of Gravity
m = 164*u.u m = DY_MASS
gravity_axis = np.array([0, 0, 1]) gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity'] tilt_gravity = options['tilt_gravity']
theta = options['theta'] theta = options['theta']
@ -489,37 +514,73 @@ def computeTrapPotential(w_x, w_z, Power, options):
R = rotation_matrix(tilt_axis, np.radians(theta)) R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis) gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis] gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = crossed_beam_potential(Positions, delta, waists, P = Power) TrappingPotential = crossed_beam_potential(Positions, waists, Power, options)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m) 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) TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif not gravity and astigmatism: elif not gravity and astigmatism:
# Influence of Astigmatism # Influence of Astigmatism
disp_foci = options['disp_foci_crossed'] TrappingPotential = astigmatic_crossed_beam_potential(Positions, waists, Power, options)
TrappingPotential = astigmatic_crossed_beam_potential(Positions, delta, waists, P = Power, del_y = disp_foci)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis])
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif gravity and astigmatism: elif gravity and astigmatism:
# Influence of Gravity and Astigmatism # Influence of Gravity and Astigmatism
m = 164*u.u m = DY_MASS
gravity_axis = np.array([0, 0, 1]) gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity'] tilt_gravity = options['tilt_gravity']
theta = options['theta'] theta = options['theta']
tilt_axis = options['tilt_axis'] tilt_axis = options['tilt_axis']
disp_foci = options['disp_foci_crossed']
if tilt_gravity: if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta)) R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis) gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis] gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = astigmatic_crossed_beam_potential(Positions, delta, waists, P = Power, del_y = disp_foci) TrappingPotential = astigmatic_crossed_beam_potential(Positions, waists, Power, options)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m) 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) TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
else: else:
TrappingPotential = IdealTrappingPotential TrappingPotential = IdealTrappingPotential
return Positions, TrappingPotential infls = np.where(np.diff(np.sign(np.gradient(np.gradient(TrappingPotential[axis].value)))))[0]
try:
if TrappingPotential[axis][0] > TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][infls[1]:-1]) - min(TrappingPotential[axis][infls[0]:infls[1]])
elif TrappingPotential[axis][0] < TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][0:infls[0]]) - min(TrappingPotential[axis][infls[0]:infls[1]])
else:
EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin
except:
EffectiveTrapDepthInKelvin = np.nan
TrapDepthsInKelvin = [IdealTrapDepthInKelvin, EffectiveTrapDepthInKelvin]
v_x = calculateCrossedBeamTrapFrequency(options['delta'], [w_x, w_z], Power, dir = 'x')
v_y = calculateCrossedBeamTrapFrequency(options['delta'], [w_x, w_z], Power, dir = 'y')
v_z = calculateCrossedBeamTrapFrequency(options['delta'], [w_x, w_z], Power, dir = 'z')
CalculatedTrapFrequencies = [v_x, v_y, v_z]
v, dv, popt, pcov = extractTrapFrequency(Positions, IdealTrappingPotential, axis)
if np.isinf(v):
v = np.nan
if np.isinf(dv):
dv = np.nan
IdealTrapFrequency = [v, dv]
if options['extract_trap_frequencies']:
v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, axis)
if np.isinf(v):
v = np.nan
if np.isinf(dv):
dv = np.nan
TrapFrequency = [v, dv]
ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency]
else:
ExtractedTrapFrequencies = [IdealTrapFrequency]
return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies
def extractWaist(Positions, TrappingPotential): def extractWaist(Positions, TrappingPotential):
tmp_pos = Positions.value tmp_pos = Positions.value
@ -538,7 +599,7 @@ def extractWaist(Positions, TrappingPotential):
popt, pcov = curve_fit(gaussian_potential, xdata, Potential, p0) popt, pcov = curve_fit(gaussian_potential, xdata, Potential, p0)
return popt, pcov return popt, pcov
def computeIntensityProfileAndPotentials(Power, waists, wavelength, options, alpha = 184.4): def computeIntensityProfileAndPotentials(Power, waists, wavelength, options, alpha = DY_POLARIZABILITY):
w_x = waists[0] w_x = waists[0]
w_z = waists[1] w_z = waists[1]
extent = options['extent'] extent = options['extent']
@ -939,4 +1000,3 @@ def plotCollisionRatesAndPSD(Gamma_elastic, PSD, modulation_depth, new_aspect_ra
##################################################################### #####################################################################
#Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability