From 4714556ab48e88983ddb048a22847908e7b23643 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 2 Mar 2023 10:08:30 +0100 Subject: [PATCH] Corrected how the effective trap depth is calculated, minor changes to code for plotting of the potentials. --- calculateDipoleTrapPotential.py | 234 +++++++++++++++++--------------- 1 file changed, 128 insertions(+), 106 deletions(-) diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py index 7d29bb7..d524b2f 100644 --- a/calculateDipoleTrapPotential.py +++ b/calculateDipoleTrapPotential.py @@ -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,11 +425,24 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): CalculatedTrapFrequencies = [v_x, v_y, v_z] v, dv, popt, pcov = extractTrapFrequency(Positions, IdealTrappingPotential, axis) - IdealTrapFrequency = [v, dv] - v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, axis) - TrapFrequency = [v, dv] - ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency] + 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 else: @@ -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] - v = Params[j][2][1][0] - dv = Params[j][2][1][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): @@ -608,7 +630,7 @@ def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterat if i % 2 == 0: plt.plot(Positions[axis], ComputedPotentials[i][axis], '--', label = 'Trap Depth = ' + str(round(IdealTrapDepthInKelvin.value, 2)) + ' ' + str(IdealTrapDepthInKelvin.unit) + '; ' + generate_label(idealv, idealdv)) elif i % 2 != 0: - plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'Effective Trap Depth = ' + str(round(EffectiveTrapDepthInKelvin.value, 2)) + ' ' + str(EffectiveTrapDepthInKelvin.unit) + '; ' + generate_label(v, dv)) + plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'Effective Trap Depth = ' + str(round(EffectiveTrapDepthInKelvin.value, 2)) + ' ' + str(EffectiveTrapDepthInKelvin.unit) + '; ' + generate_label(v, dv)) else: if i % 2 == 0: plt.plot(Positions[axis], ComputedPotentials[i][axis], '--', label = 'Trap Depth = ' + str(round(IdealTrapDepthInKelvin.value, 2)) + ' ' + str(IdealTrapDepthInKelvin.unit) + '; ' + generate_label(idealv, idealdv)) @@ -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 = [] + + # """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]) + 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, 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 = { @@ -1058,85 +1082,85 @@ if __name__ == '__main__': # 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))) + # 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 - 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) + # 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))) - v_x = np.zeros(len(w_x)) - #v_y = np.zeros(len(w_x)) - v_z = np.zeros(len(w_x)) + # 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) - 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 = np.zeros(len(w_x)) + # #v_y = np.zeros(len(w_x)) + # v_z = np.zeros(len(w_x)) - #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 + # 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 - 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] + # #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 - calc_alpha_deviation = [(g - h) for g, h in zip(alpha_x, alpha_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] - 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() + # calc_alpha_deviation = [(g - h) for g, h in zip(alpha_x, alpha_z)] - 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, 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 @@ -1180,6 +1204,4 @@ if __name__ == '__main__': # plt.xlim([-500, 500]) # plt.ylim([-1800, -200]) # plt.legend() - # plt.show() - - \ No newline at end of file + # plt.show() \ No newline at end of file