MAJOR rewrite - rewrote code snippets that plotted from measurement data, added calculations of trap frequencies for crossed dipole trap config, put running of code in to a Jupyter notebook.

This commit is contained in:
Karthik 2023-03-14 17:13:37 +01:00
parent 4714556ab4
commit 2fdf592dda
2 changed files with 1070 additions and 391 deletions

1006
ODTCalculations.ipynb Normal file

File diff suppressed because one or more lines are too long

View File

@ -80,64 +80,43 @@ 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):
sigma = 8.468e-14 * (u.m)**(2)
return (1/eta * 1/sigma * NCount * pixel_size**2/magnification**2).decompose()
def meanThermalVelocity(T, m = 164*u.u): def meanThermalVelocity(T, m = 164*u.u):
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, Polarizability, N, T, m = 164*u.u, use_measured_tf = False): # For a thermal cloud def particleDensity(w_x, w_z, Power, Polarizability, N, T, m = 164*u.u): # For a thermal cloud
if not use_measured_tf: v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x')
v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x') v_y = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'y')
v_y = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'y') v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'z')
v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, 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)
else:
fin_mod_dep = [0.5, 0.3, 0.7, 0.9, 0.8, 1.0, 0.6, 0.4, 0.2, 0.1]
v_x = [0.28, 0.690, 0.152, 0.102, 0.127, 0.099, 0.205, 0.404, 1.441, 2.813] * u.kHz
dv_x = [0.006, 0.005, 0.006, 0.003, 0.002, 0.002,0.002, 0.003, 0.006, 0.024] * u.kHz
v_z = [1.278, 1.719, 1.058, 0.923, 0.994, 0.911, 1.157, 1.446, 2.191, 2.643] * u.kHz
dv_z = [0.007, 0.009, 0.007, 0.005, 0.004, 0.004, 0.005, 0.007, 0.009, 0.033] * u.kHz
sorted_fin_mod_dep, sorted_v_x = zip(*sorted(zip(fin_mod_dep, v_x)))
sorted_fin_mod_dep, sorted_dv_x = zip(*sorted(zip(fin_mod_dep, dv_x)))
sorted_fin_mod_dep, sorted_v_z = zip(*sorted(zip(fin_mod_dep, v_z)))
sorted_fin_mod_dep, sorted_dv_z = zip(*sorted(zip(fin_mod_dep, dv_z)))
fin_mod_dep = [1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]
v_y = [3.08, 3.13, 3.27, 3.46, 3.61, 3.82, 3.51, 3.15, 3.11, 3.02] * u.Hz
dv_y = [0.03, 0.04, 0.04, 0.05, 0.07, 0.06, 0.11, 0.07, 0.1, 1.31] * u.Hz
sorted_fin_mod_dep, sorted_v_y = zip(*sorted(zip(fin_mod_dep, v_y)))
sorted_fin_mod_dep, sorted_dv_y = zip(*sorted(zip(fin_mod_dep, dv_y)))
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):
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]
dalpha_y = [alpha_y[i] * np.sqrt((dv_z[0]/v_z[0])**2 + (dv_z[i]/v_z[i])**2) for i in range(len(v_z))] dalpha_y = [alpha_y[i] * np.sqrt((dv_z[0]/v_z[0])**2 + (dv_z[i]/v_z[i])**2) for i in range(len(v_z))]
avg_alpha = [(g + h) / 2 for g, h in zip(alpha_x, alpha_y)] avg_alpha = [(g + h) / 2 for g, h in zip(alpha_x, alpha_y)]
sorted_fin_mod_dep, new_aspect_ratio = zip(*sorted(zip(fin_mod_dep, (w_x * avg_alpha) / w_z))) new_aspect_ratio = (w_x * avg_alpha) / w_z
fin_mod_dep = [1.0, 0.8, 0.6, 0.4, 0.2, 0.9, 0.7, 0.5, 0.3, 0.1]
T_x = [22.1, 27.9, 31.7, 42.2, 145.8, 27.9, 33.8, 42.4, 61.9, 136.1] * u.uK
dT_x = [1.7, 2.6, 2.4, 3.7, 1.1, 2.2, 3.2, 1.7, 2.2, 1.2] * u.uK
T_y = [13.13, 14.75, 18.44, 26.31, 52.55, 13.54, 16.11, 21.15, 35.81, 85.8] * u.uK
dT_y = [0.05, 0.05, 0.07, 0.16, 0.28, 0.04, 0.07, 0.10, 0.21, 0.8] * u.uK
avg_T = [(g + h) / 2 for g, h in zip(T_x, T_y)] avg_T = [(g + h) / 2 for g, h in zip(T_x, T_y)]
avg_dT = [0.5 * np.sqrt(g**2 + h**2) for g, h in zip(dT_x, dT_y)] avg_dT = [0.5 * np.sqrt(g**2 + h**2) for g, h in zip(dT_x, dT_y)]
sorted_fin_mod_dep, sorted_avg_T = zip(*sorted(zip(fin_mod_dep, avg_T))) pd = np.zeros(len(modulation_depth))
sorted_fin_mod_dep, sorted_avg_dT = zip(*sorted(zip(fin_mod_dep, avg_dT))) dpd = np.zeros(len(modulation_depth))
pd = np.zeros(len(fin_mod_dep))
dpd = np.zeros(len(fin_mod_dep))
for i in range(len(fin_mod_dep)): for i in range(len(modulation_depth)):
particle_density = (N * (2 * np.pi)**3 * (sorted_v_x[i] * sorted_v_y[i] * sorted_v_z[i]) * (m / (2 * np.pi * ac.k_B * sorted_avg_T[i]))**(3/2)).decompose() particle_density = (N * (2 * np.pi)**3 * (v_x[i] * v_y[i] * v_z[i]) * (m / (2 * np.pi * ac.k_B * avg_T[i]))**(3/2)).decompose()
pd[i] = particle_density.value pd[i] = particle_density.value
dpd[i] = (((N * (2 * np.pi)**3 * (m / (2 * np.pi * ac.k_B * sorted_avg_T[i]))**(3/2)) * ((sorted_dv_x[i] * sorted_v_y[i] * sorted_v_z[i]) + (sorted_v_x[i] * sorted_dv_y[i] * sorted_v_z[i]) + (sorted_v_x[i] * sorted_v_y[i] * sorted_dv_z[i]) - (1.5*(sorted_v_x[i] * sorted_v_y[i] * sorted_v_z[i])*(sorted_avg_dT[i]/sorted_avg_T[i])))).decompose()).value dpd[i] = (((N * (2 * np.pi)**3 * (m / (2 * np.pi * ac.k_B * avg_T[i]))**(3/2)) * ((dv_x[i] * v_y[i] * v_z[i]) + (v_x[i] * dv_y[i] * v_z[i]) + (v_x[i] * v_y[i] * dv_z[i]) - (1.5*(v_x[i] * v_y[i] * v_z[i])*(avg_dT[i]/avg_T[i])))).decompose()).value
pd = pd*particle_density.unit pd = pd*particle_density.unit
dpd = dpd*particle_density.unit dpd = dpd*particle_density.unit
return pd, dpd, sorted_avg_T, sorted_avg_dT, new_aspect_ratio, sorted_fin_mod_dep return pd, dpd, avg_T, avg_dT, new_aspect_ratio
def thermaldeBroglieWavelength(T, m = 164*u.u): def thermaldeBroglieWavelength(T, m = 164*u.u):
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))
@ -290,7 +269,7 @@ def crossed_beam_potential(positions, theta, waists, P, alpha, wavelength=1.064*
def trap_depth(w_1, w_2, P, alpha): def trap_depth(w_1, w_2, P, alpha):
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, Polarizability, m = 164*u.u, dir = 'x'): def calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x', m = 164*u.u):
TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability) TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability)
TrapFrequency = np.nan TrapFrequency = np.nan
if dir == 'x': if dir == 'x':
@ -302,15 +281,48 @@ def calculateTrapFrequency(w_x, w_z, Power, Polarizability, m = 164*u.u, dir = '
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, Polarizability = 184.4, wavelength=1.064*u.um):
TrapDepth_1 = trap_depth(Waists[0][0], Waists[1][0], Powers[0], alpha=Polarizability)
TrapDepth_2 = trap_depth(Waists[0][1], Waists[1][1], Powers[1], alpha=Polarizability)
w_x1 = Waists[0][0]
w_z1 = Waists[1][0]
w_y2 = Waists[0][1]
w_z2 = Waists[1][1]
zReff_1 = np.sqrt(2) * z_R(w_x1, 1.064*u.um) * z_R(w_z1, 1.064*u.um) / np.sqrt(z_R(w_x1, 1.064*u.um)**2 + z_R(w_z1, 1.064*u.um)**2)
zReff_2 = np.sqrt(2) * z_R(w_y2, 1.064*u.um) * z_R(w_z2, 1.064*u.um) / np.sqrt(z_R(w_y2, 1.064*u.um)**2 + z_R(w_z2, 1.064*u.um)**2)
wy2alpha = np.sqrt(((np.cos(np.radians(90 - delta)) / w_y2)**2 + (np.sin(np.radians(90 - delta)) / (2 * zReff_2))**2)**(-1))
zR2alpha = np.sqrt(((np.sin(np.radians(90 - delta)) / w_y2)**2 + (np.cos(np.radians(90 - delta)) / (2 * zReff_2))**2)**(-1))
TrapFrequency = np.nan
if dir == 'x':
v_1 = (1/(2 * np.pi)) * np.sqrt(4/m * (TrapDepth_1 / w_x1**2))
v_2 = (1/(2 * np.pi)) * np.sqrt(4/m * (TrapDepth_2 / zR2alpha**2))
TrapFrequency = (np.sqrt(v_1**2 + v_2**2)).decompose()
elif dir == 'y':
v_1 = (1/(2 * np.pi)) * np.sqrt(1/m * (2 * TrapDepth_1/ zReff_1**2))
v_2 = (1/(2 * np.pi)) * np.sqrt(1/m * (4 * TrapDepth_2/ wy2alpha**2))
TrapFrequency = (np.sqrt(v_1**2 + v_2**2)).decompose()
elif dir == 'z':
v_1 = (1/(2 * np.pi)) * np.sqrt(4/m * (TrapDepth_1 / w_z1**2))
v_2 = (1/(2 * np.pi)) * np.sqrt(4/m * (TrapDepth_2 / w_z2**2))
TrapFrequency = (np.sqrt(v_1**2 + v_2**2)).decompose()
return round(TrapFrequency.value, 2)*u.Hz
def extractTrapFrequency(Positions, TrappingPotential, axis): def extractTrapFrequency(Positions, TrappingPotential, axis):
tmp_pos = Positions[axis, :] tmp_pos = Positions[axis, :]
tmp_pot = TrappingPotential[axis] tmp_pot = TrappingPotential[axis]
center_idx = np.argmin(tmp_pot) center_idx = np.argmin(tmp_pot)
lb = int(round(center_idx - len(tmp_pot)/500, 1)) lb = int(round(center_idx - len(tmp_pot)/200, 1))
ub = int(round(center_idx + len(tmp_pot)/500, 1)) ub = int(round(center_idx + len(tmp_pot)/200, 1))
xdata = tmp_pos[lb:ub] xdata = tmp_pos[lb:ub]
Potential = tmp_pot[lb:ub] Potential = tmp_pot[lb:ub]
p0=[1e3, tmp_pos[center_idx].value, -100] p0 = [1e3, tmp_pos[center_idx].value, -100]
popt, pcov = curve_fit(harmonic_potential, xdata, Potential, p0) popt, pcov = curve_fit(harmonic_potential, xdata, Potential, p0)
v = popt[0] v = popt[0]
dv = pcov[0][0]**0.5 dv = pcov[0][0]**0.5
@ -359,9 +371,9 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
IdealTrappingPotential = (IdealTrappingPotential/ac.k_B).to(u.uK) IdealTrappingPotential = (IdealTrappingPotential/ac.k_B).to(u.uK)
else: else:
theta = options['delta'] 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, theta, waists, P = Power, alpha = Polarizability) IdealTrappingPotential = crossed_beam_potential(Positions, delta, waists, P = Power, alpha = Polarizability)
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)
@ -758,17 +770,8 @@ def plotTrapFrequencies(v_x, v_y, v_z, modulation_depth, new_aspect_ratio, plot_
ax3.legend(lns, labs, prop={'size': 12, 'weight': 'bold'}) ax3.legend(lns, labs, prop={'size': 12, 'weight': 'bold'})
plt.show() plt.show()
def plotMeasuredTrapFrequencies(w_x, w_z, plot_against_mod_depth = True): def plotMeasuredTrapFrequencies(fx, dfx, fy, dfy, fz, dfz, modulation_depth_radial, modulation_depth_axial, w_x, w_z, plot_against_mod_depth = True):
fin_mod_dep = [0, 0.5, 0.3, 0.7, 0.9, 0.8, 1.0, 0.6, 0.4, 0.2, 0.1]
fx = [3.135, 0.28, 0.690, 0.152, 0.102, 0.127, 0.099, 0.205, 0.404, 1.441, 2.813]
dfx = [0.016, 0.006, 0.005, 0.006, 0.003, 0.002, 0.002,0.002, 0.003, 0.006, 0.024]
fz = [2.746, 1.278, 1.719, 1.058, 0.923, 0.994, 0.911, 1.157, 1.446, 2.191, 2.643]
dfz = [0.014, 0.007, 0.009, 0.007, 0.005, 0.004, 0.004, 0.005, 0.007, 0.009, 0.033]
fin_mod_dep_y = [1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]
fy = [3.08, 3.13, 3.27, 3.46, 3.61, 3.82, 3.51, 3.15, 3.11, 3.02]
dfy = [0.03, 0.04, 0.04, 0.05, 0.07, 0.06, 0.11, 0.07, 0.1, 1.31]
alpha_x = [(fx[0]/x)**(2/3) for x in fx] alpha_x = [(fx[0]/x)**(2/3) for x in fx]
dalpha_x = [alpha_x[i] * np.sqrt((dfx[0]/fx[0])**2 + (dfx[i]/fx[i])**2) for i in range(len(fx))] dalpha_x = [alpha_x[i] * np.sqrt((dfx[0]/fx[0])**2 + (dfx[i]/fx[i])**2) for i in range(len(fx))]
alpha_y = [(fy[0]/y)**2 for y in fy] alpha_y = [(fy[0]/y)**2 for y in fy]
@ -781,9 +784,9 @@ def plotMeasuredTrapFrequencies(w_x, w_z, plot_against_mod_depth = True):
if plot_against_mod_depth: if plot_against_mod_depth:
fig, ax1 = plt.subplots(figsize=(8, 6)) fig, ax1 = plt.subplots(figsize=(8, 6))
ax2 = ax1.twinx() ax2 = ax1.twinx()
ax1.errorbar(fin_mod_dep, fx, yerr = dfx, fmt= 'or', label = 'v_x', markersize=5, capsize=5) ax1.errorbar(modulation_depth_radial, fx, yerr = dfx, fmt= 'or', label = 'v_x', markersize=5, capsize=5)
ax2.errorbar(fin_mod_dep_y, fy, yerr = dfy, fmt= '*g', label = 'v_y', markersize=5, capsize=5) ax2.errorbar(modulation_depth_axial, fy, yerr = dfy, fmt= '*g', label = 'v_y', markersize=5, capsize=5)
ax1.errorbar(fin_mod_dep, fz, yerr = dfz, fmt= '^b', label = 'v_z', markersize=5, capsize=5) ax1.errorbar(modulation_depth_radial, fz, yerr = dfz, fmt= '^b', label = 'v_z', markersize=5, capsize=5)
ax1.set_xlabel('Modulation depth', fontsize= 12, fontweight='bold') ax1.set_xlabel('Modulation depth', fontsize= 12, fontweight='bold')
ax1.set_ylabel('Trap Frequency (kHz)', fontsize= 12, fontweight='bold') ax1.set_ylabel('Trap Frequency (kHz)', fontsize= 12, fontweight='bold')
ax1.tick_params(axis="y", labelcolor='b') ax1.tick_params(axis="y", labelcolor='b')
@ -804,27 +807,11 @@ def plotMeasuredTrapFrequencies(w_x, w_z, plot_against_mod_depth = True):
plt.grid(visible=1) plt.grid(visible=1)
plt.show() plt.show()
def plotRatioOfTrapFrequencies(plot_against_mod_depth = True): def plotRatioOfTrapFrequencies(fx, fy, fz, dfx, dfy, dfz, v_x, v_y, v_z, modulation_depth, w_x, w_z, plot_against_mod_depth = True):
modulation_depth = [0.5, 0.3, 0.7, 0.9, 0.8, 1.0, 0.6, 0.4, 0.2, 0.1]
w_xs = w_x * convert_modulation_depth_to_alpha(modulation_depth)[0] w_xs = w_x * convert_modulation_depth_to_alpha(modulation_depth)[0]
new_aspect_ratio = w_xs / w_z new_aspect_ratio = w_xs / w_z
v_x = np.zeros(len(modulation_depth))
v_y = np.zeros(len(modulation_depth))
v_z = np.zeros(len(modulation_depth))
for i in range(len(modulation_depth)):
v_x[i] = calculateTrapFrequency(w_xs[i], w_z, Power, Polarizability, dir = 'x').value / 1e3
v_y[i] = calculateTrapFrequency(w_xs[i], w_z, Power, Polarizability, dir = 'y').value
v_z[i] = calculateTrapFrequency(w_xs[i], w_z, Power, Polarizability, dir = 'z').value / 1e3
fx = [0.28, 0.690, 0.152, 0.102, 0.127, 0.099, 0.205, 0.404, 1.441, 2.813]
dfx = [0.006, 0.005, 0.006, 0.003, 0.002, 0.002,0.002, 0.003, 0.006, 0.024]
fy = [3.08, 3.13, 3.27, 3.46, 3.61, 3.82, 3.51, 3.15, 3.11, 3.02]
dfy = [0.03, 0.04, 0.04, 0.05, 0.07, 0.06, 0.11, 0.07, 0.1, 1.31]
fz = [1.278, 1.719, 1.058, 0.923, 0.994, 0.911, 1.157, 1.446, 2.191, 2.643]
dfz = [0.007, 0.009, 0.007, 0.005, 0.004, 0.004, 0.005, 0.007, 0.009, 0.033]
plt.figure() plt.figure()
if plot_against_mod_depth: if plot_against_mod_depth:
@ -890,318 +877,4 @@ def plotCollisionRatesAndPSD(Gamma_elastic, PSD, modulation_depth, new_aspect_ra
plt.grid(visible=1) plt.grid(visible=1)
plt.show() plt.show()
##################################################################### #####################################################################
# RUN SCRIPT WITH OPTIONS BELOW #
#####################################################################
if __name__ == '__main__':
Power = 0.420*u.W
Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability
Wavelength = 1.064*u.um
w_x, w_z = 30*u.um, 30*u.um # Beam Waists in the x and y directions
# Power = 11*u.W
# Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability
# w_x, w_z = 54.0*u.um, 54.0*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
'crossed': False,
'delta': 70, # angle between arms in degrees
'modulation': False,
'aspect_ratio': 4, # required aspect ratio of modulated arm
'gravity': True,
'tilt_gravity': True,
'theta': 0.1, # gravity tilt angle in degrees
'tilt_axis': [1, 0, 0], # lab space coordinates are rotated about x-axis in reference frame of beam
'astigmatism': False,
'disp_foci': 0.9 * z_R(w_0 = np.asarray([30]), lamb = 1.064)[0]*u.um, # difference in position of the foci along the propagation direction (Astigmatism)
'extract_trap_frequencies': True
}
# """Plot ideal trap potential resulting for given parameters only"""
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])
ComputedPotentials = np.asarray(ComputedPotentials)
plotPotential(Positions, ComputedPotentials, options, Params)
"""Plot harmonic fit for trap potential resulting for given parameters only"""
# v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, options['axis'])
# plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, options['axis'], popt, pcov)
"""Plot trap potential resulting for given parameters (with one parameter being a list of values and the potential being computed for each of these values) only"""
# ComputedPotentials = []
# Params = []
# Power = [10, 20, 25, 30, 35, 40]*u.W # Single Beam Power
# for p in Power:
# Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies = computeTrapPotential(w_x, w_z, p, Polarizability, options)
# ComputedPotentials.append(IdealTrappingPotential)
# ComputedPotentials.append(TrappingPotential)
# Params.append([TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies])
# ComputedPotentials = np.asarray(ComputedPotentials)
# plotPotential(Positions, ComputedPotentials, options, Params)
"""Plot transverse intensity profile and trap potential resulting for given parameters only"""
# options = {
# 'extent': 60, # range of spatial coordinates in one direction to calculate trap potential over
# 'modulation': True,
# 'modulation_function': 'arccos',
# 'modulation_amplitude': 2.16
# }
# positions, waists, I, U, p = computeIntensityProfileAndPotentials(Power, [w_x, w_z], Polarizability, Wavelength, options)
# plotIntensityProfileAndPotentials(positions, waists, I, U)
"""Plot gaussian fit for trap potential resulting from modulation for given parameters only"""
# x_Positions = positions[0].value
# z_Positions = positions[1].value
# x_Potential = U[:, np.where(z_Positions==0)[0][0]].value
# z_Potential = U[np.where(x_Positions==0)[0][0], :].value
# poptx, pcovx = p[0], p[1]
# poptz, pcovz = p[2], p[3]
# plotGaussianFit(x_Positions, x_Potential, poptx, pcovx)
# plotGaussianFit(z_Positions, z_Potential, poptz, pcovz)
"""Calculate relevant parameters for evaporative cooling"""
# AtomNumber = 1.00 * 1e7
# BField = 2.5 * u.G
# modulation = True
# if modulation:
# modulation_depth = 0.6
# w_x = w_x * convert_modulation_depth_to_alpha(modulation_depth)[0]
# Temperature = convert_modulation_depth_to_temperature(modulation_depth)[0] * u.uK
# else:
# modulation_depth = 0.0
# Temperature = convert_modulation_depth_to_temperature(modulation_depth)[0] * u.uK
# 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('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 = %.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))
# print('a_s = %.2f ' %(scatteringLength(BField)[0] / ac.a0))
"""Calculate relevant parameters for evaporative cooling for different modulation depths, temperatures"""
# AtomNumber = 1.00 * 1e7
# BField = 1.4 * u.G
# modulation_depth = np.arange(0, 1.0, 0.08)
# w_xs = w_x * convert_modulation_depth_to_alpha(modulation_depth)[0]
# new_aspect_ratio = w_xs / w_z
# Temperatures = convert_modulation_depth_to_temperature(modulation_depth)[0] * u.uK
# plot_against_mod_depth = True
# # n = np.zeros(len(modulation_depth))
# Gamma_elastic = np.zeros(len(modulation_depth))
# PSD = np.zeros(len(modulation_depth))
# v_x = np.zeros(len(modulation_depth))
# v_y = np.zeros(len(modulation_depth))
# v_z = np.zeros(len(modulation_depth))
# for i in range(len(modulation_depth)):
# # n[i] = particleDensity(w_xs[i], w_z, Power, Polarizability, N = AtomNumber, T = Temperatures[i], m = 164*u.u).decompose().to(u.cm**(-3))
# Gamma_elastic[i] = calculateElasticCollisionRate(w_xs[i], w_z, Power, Polarizability, N = AtomNumber, T = Temperatures[i], B = BField).value
# PSD[i] = calculatePSD(w_xs[i], w_z, Power, Polarizability, N = AtomNumber, T = Temperatures[i]).decompose().value
# v_x[i] = calculateTrapFrequency(w_xs[i], w_z, Power, Polarizability, dir = 'x').value
# v_y[i] = calculateTrapFrequency(w_xs[i], w_z, Power, Polarizability, dir = 'y').value
# v_z[i] = calculateTrapFrequency(w_xs[i], w_z, Power, Polarizability, dir = 'z').value
"""Plot alphas"""
# plotAlphas()
"""Plot Temperatures"""
# plotTemperatures(w_x, w_z, plot_against_mod_depth = plot_against_mod_depth)
"""Plot trap frequencies"""
# plotTrapFrequencies(v_x, v_y, v_z, modulation_depth, new_aspect_ratio, plot_against_mod_depth = plot_against_mod_depth)
# plotMeasuredTrapFrequencies(w_x, w_z, plot_against_mod_depth = plot_against_mod_depth)
# plotRatioOfTrapFrequencies(plot_against_mod_depth = plot_against_mod_depth)
"""Plot Feshbach Resonances"""
# plotScatteringLengths(B_range = [0, 3.6])
"""Plot Collision Rates and PSD"""
# plotCollisionRatesAndPSD(Gamma_elastic, PSD, modulation_depth, new_aspect_ratio, plot_against_mod_depth = plot_against_mod_depth)
"""Plot Collision Rates and PSD from only measured trap frequencies"""
# pd, dpd, T, dT, new_aspect_ratio, modulation_depth = particleDensity(w_x, w_z, Power, Polarizability, AtomNumber, 0, m = 164*u.u, use_measured_tf = True)
# Gamma_elastic = [(pd[i] * scatteringCrossSection(BField) * meanThermalVelocity(T[i]) / (2 * np.sqrt(2))).decompose() for i in range(len(pd))]
# Gamma_elastic_values = [(Gamma_elastic[i]).value for i in range(len(Gamma_elastic))]
# dGamma_elastic = [(Gamma_elastic[i] * ((dpd[i]/pd[i]) + (dT[i]/(2*T[i])))).decompose() for i in range(len(Gamma_elastic))]
# dGamma_elastic_values = [(dGamma_elastic[i]).value for i in range(len(dGamma_elastic))]
# PSD = [((pd[i] * thermaldeBroglieWavelength(T[i])**3).decompose()).value for i in range(len(pd))]
# dPSD = [((PSD[i] * ((dpd[i]/pd[i]) - (1.5 * dT[i]/T[i]))).decompose()).value for i in range(len(Gamma_elastic))]
# fig, ax1 = plt.subplots(figsize=(8, 6))
# ax2 = ax1.twinx()
# ax1.errorbar(modulation_depth, Gamma_elastic_values, yerr = dGamma_elastic_values, fmt = 'ob', markersize=5, capsize=5)
# ax2.errorbar(modulation_depth, PSD, yerr = dPSD, fmt = '-^r', markersize=5, capsize=5)
# ax2.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
# ax1.set_xlabel('Modulation depth', fontsize= 12, fontweight='bold')
# ax1.set_ylabel('Elastic Collision Rate (' + str(Gamma_elastic[0].unit) + ')', fontsize= 12, fontweight='bold')
# ax1.tick_params(axis="y", labelcolor='b')
# ax2.set_ylabel('Phase Space Density', fontsize= 12, fontweight='bold')
# ax2.tick_params(axis="y", labelcolor='r')
# plt.tight_layout()
# plt.grid(visible=1)
# plt.show()
""" Investigate deviation in alpha"""
# Power = 40*u.W
# Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability
# Wavelength = 1.064*u.um
# w_x, w_z = 27.5*u.um, 33.8*u.um
# options = {
# '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
# 'crossed': False,
# 'theta': 0,
# 'modulation': False,
# 'gravity': True,
# 'tilt_gravity': True,
# 'theta': 10, # in degrees
# 'tilt_axis': [1, 0, 0], # lab space coordinates are rotated about x-axis in reference frame of beam
# 'astigmatism': True,
# 'disp_foci': 0.9 * z_R(w_0 = np.asarray([30]), lamb = 1.064)[0]*u.um # difference in position of the foci along the propagation direction (Astigmatism)
# }
# modulation_depth = np.arange(0, 1.1, 0.1)
# Alphas, fin_mod_dep, meas_alpha_x, meas_alpha_z, dalpha_x, dalpha_z = convert_modulation_depth_to_alpha(modulation_depth)
# meas_alpha_deviation = [(g - h) for g, h in zip(meas_alpha_x, meas_alpha_z)]
# sorted_fin_mod_dep, sorted_meas_alpha_deviation = zip(*sorted(zip(fin_mod_dep, meas_alpha_deviation)))
# avg_alpha = [(g + h) / 2 for g, h in zip(meas_alpha_x, meas_alpha_z)]
# sorted_fin_mod_dep, new_aspect_ratio = zip(*sorted(zip(fin_mod_dep, (w_x * avg_alpha) / w_z)))
# current_ar = w_x/w_z
# aspect_ratio = np.arange(current_ar, 10*current_ar, 0.8)
# w_x = w_x * (aspect_ratio / current_ar)
# v_x = np.zeros(len(w_x))
# #v_y = np.zeros(len(w_x))
# v_z = np.zeros(len(w_x))
# for i in range(len(w_x)):
# options['axis'] = 0
# ExtractedTrapFrequencies = computeTrapPotential(w_x[i], w_z, Power, Polarizability, options)[5]
# tmp = ExtractedTrapFrequencies[1][0]
# v_x[i] = tmp if not np.isinf(tmp) else np.nan
# #options['axis'] = 1
# #ExtractedTrapFrequencies = computeTrapPotential(w_x[i], w_z, Power, Polarizability, options)[5]
# #tmp = ExtractedTrapFrequencies[1][0]
# #v_y[i] = tmp if not np.isinf(tmp) else np.nan
# options['axis'] = 2
# ExtractedTrapFrequencies = computeTrapPotential(w_x[i], w_z, Power, Polarizability, options)[5]
# tmp = ExtractedTrapFrequencies[1][0]
# v_z[i] = tmp if not np.isinf(tmp) else np.nan
# #v_x[i] = calculateTrapFrequency(w_x[i], w_z, Power, Polarizability, dir = 'x').value
# #v_y[i] = calculateTrapFrequency(w_x[i], w_z, Power, Polarizability, dir = 'y').value
# #v_z[i] = calculateTrapFrequency(w_x[i], w_z, Power, Polarizability, dir = 'z').value
# alpha_x = [(v_x[0]/v)**(2/3) for v in v_x]
# alpha_z = [(v_z[0]/v)**2 for v in v_z]
# calc_alpha_deviation = [(g - h) for g, h in zip(alpha_x, alpha_z)]
# plt.figure()
# plt.plot(aspect_ratio, alpha_x, '-o', label = 'From horz TF')
# plt.plot(aspect_ratio, alpha_z, '-^', label = 'From vert TF')
# plt.xlabel('Aspect Ratio', fontsize= 12, fontweight='bold')
# plt.ylabel('$\\alpha$', fontsize= 12, fontweight='bold')
# plt.tight_layout()
# plt.grid(visible=1)
# plt.legend(prop={'size': 12, 'weight': 'bold'})
# plt.show()
# plt.figure()
# plt.plot(aspect_ratio, calc_alpha_deviation, '--ob', label = 'Extracted')
# plt.plot(new_aspect_ratio, sorted_meas_alpha_deviation, '-or', label = 'Measured')
# plt.xlabel('Aspect Ratio', fontsize= 12, fontweight='bold')
# plt.ylabel('$\\Delta \\alpha$', fontsize= 12, fontweight='bold')
# plt.ylim([-0.5, 1])
# plt.tight_layout()
# plt.grid(visible=1)
# plt.legend(prop={'size': 12, 'weight': 'bold'})
# plt.show()
"""Plot ideal crossed beam trap potential resulting for given parameters only"""
# Powers = [40, 10] * u.W
# Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability
# Wavelength = 1.064*u.um
# w_x = [27.5, 54]*u.um # Beam Waists in the x direction
# w_z = [33.8, 54]*u.um # Beam Waists in the y direction
# Powers = [30, 8] * u.W
# Polarizability = 136 # in a.u, most precise measured value of Dy polarizability
# Wavelength = 1.064*u.um
# w_x = [20.5, 101.3]*u.um # Beam Waists in the x direction
# w_z = [20.5, 95.0]*u.um # Beam Waists in the y direction
# options = {
# 'axis': 1, # axis referenced to the beam along which you want the dipole trap potential
# 'extent': 2e3, # range of spatial coordinates in one direction to calculate trap potential over
# 'crossed': True,
# 'delta': 70,
# 'modulation': False,
# 'aspect_ratio': 5,
# '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,
# '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)
# }
# Positions, TrapPotential = computeTrapPotential(w_x, w_z, Powers, Polarizability, options)
# plt.rcParams["figure.figsize"] = [7.00, 3.50]
# plt.rcParams["figure.autolayout"] = True
# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
# ax.scatter(TrapPotential[0], TrapPotential[1], TrapPotential[2], c=TrapPotential[2], alpha=1)
# plt.show()
# plt.figure()
# plt.plot(Positions[options['axis']], TrapPotential[options['axis']], label = 'Crossed beam potential')
# plt.xlim([-500, 500])
# plt.ylim([-1800, -200])
# plt.legend()
# plt.show()