Corrected how the effective trap depth is calculated, minor changes to code for plotting of the potentials.

This commit is contained in:
Karthik 2023-03-02 10:08:30 +01:00
parent 69ba21b1e6
commit 4714556ab4

View File

@ -407,10 +407,13 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
TrappingPotential = IdealTrappingPotential
if not crossed:
infls = np.where(np.diff(np.sign(np.gradient(np.gradient(TrappingPotential[axis].value)))))[0]
if TrappingPotential[axis][0] > TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = TrappingPotential[axis][-1] - min(TrappingPotential[axis])
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][infls[1]:-1]) - min(TrappingPotential[axis][infls[0]:infls[1]])
elif TrappingPotential[axis][0] < TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = TrappingPotential[axis][0] - min(TrappingPotential[axis])
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][0:infls[0]]) - min(TrappingPotential[axis][infls[0]:infls[1]])
else:
EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin
@ -422,10 +425,23 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
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
@ -582,7 +598,9 @@ def plotGaussianFit(Positions, TrappingPotential, popt, pcov):
plt.tight_layout()
plt.show()
def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterateOver = [], save = False):
def plotPotential(Positions, ComputedPotentials, options, Params = [], listToIterateOver = [], save = False):
axis = options['axis']
plt.figure(figsize=(9, 7))
for i in range(np.size(ComputedPotentials, 0)):
@ -598,8 +616,12 @@ def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterat
idealv = Params[j][2][0][0]
idealdv = Params[j][2][0][1]
if options['extract_trap_frequencies']:
v = Params[j][2][1][0]
dv = Params[j][2][1][1]
else:
v = np.nan
dv = np.nan
if listToIterateOver:
if np.size(ComputedPotentials, 0) == len(listToIterateOver):
@ -874,41 +896,43 @@ def plotCollisionRatesAndPSD(Gamma_elastic, PSD, modulation_depth, new_aspect_ra
if __name__ == '__main__':
# 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 # Beam Waists in the x and y directions
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': 2, # 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
'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,
'theta': 0,
'modulation': True,
'aspect_ratio': 1,
'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.5, # in degrees
'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': 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)
'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])
# """Plot ideal trap potential resulting for given parameters only"""
# ComputedPotentials = np.asarray(ComputedPotentials)
# plotPotential(Positions, ComputedPotentials, options['axis'], Params)
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'])
@ -925,7 +949,7 @@ if __name__ == '__main__':
# Params.append([TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies])
# ComputedPotentials = np.asarray(ComputedPotentials)
# plotPotential(Positions, ComputedPotentials, options['axis'], Params)
# plotPotential(Positions, ComputedPotentials, options, Params)
"""Plot transverse intensity profile and trap potential resulting for given parameters only"""
# options = {
@ -1059,83 +1083,83 @@ if __name__ == '__main__':
""" 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
# 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)
}
# 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)))
# 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)
# 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))
# 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
# 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_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] = 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
# #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]
# 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)]
# 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, 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()
# 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"""
@ -1181,5 +1205,3 @@ if __name__ == '__main__':
# plt.ylim([-1800, -200])
# plt.legend()
# plt.show()