From 7a6c4c82e4531254aecb1b7c9f1273418b39f827 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 8 Feb 2023 15:58:44 +0100 Subject: [PATCH] MAJOR CHANGE: Refactored code, added functions to calculate PSD, collision rate --- calculateDipoleTrapPotential.py | 330 ++++++++++++++++++++------------ 1 file changed, 212 insertions(+), 118 deletions(-) diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py index 6600848..f7fc8b8 100644 --- a/calculateDipoleTrapPotential.py +++ b/calculateDipoleTrapPotential.py @@ -66,6 +66,9 @@ def calculateTrapFrequency(w_x, w_z, Power, Polarizability, m = 164*u.u, dir = ' TrapFrequency = np.nan if dir == 'x': TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth / (m*w_x**2))).decompose() + elif dir == 'y': + zReff = np.sqrt(2) * z_R(w_x, 1.064*u.um) * z_R(w_z, 1.064*u.um) / np.sqrt(z_R(w_x, 1.064*u.um)**2 + z_R(w_z, 1.064*u.um)**2) + TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(2 * TrapDepth/ (m*zReff**2))).decompose() elif dir == 'z': TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth/ (m*w_z**2))).decompose() return round(TrapFrequency.value, 2)*u.Hz @@ -84,6 +87,134 @@ def extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, axis): dv = pcov[0][0]**0.5 return v, dv, popt, pcov +def meanThermalVelocity(T, m = 164*u.u): + 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): # For a thermal cloud + 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') + + return N * (2 * np.pi)**3 * (v_x * v_y * v_z) * (m / (2 * np.pi * ac.k_B * T))**(3/2) + +def thermaldeBroglieWavelength(T, m = 164*u.u): + return np.sqrt((2*np.pi*ac.hbar**2)/(m*ac.k_B*T)) + +def scatteringLength(B): + a_bkg = 87 * ac.a0 + #resonanceWidth = 0.005 * u.G + #resonanceB = 0.5 * u.G + + #return a_bkg * (1 - resonanceWidth/(B - resonanceB)) + return a_bkg + +def dipolarLength(mu = 9.93 * ac.muB, m = 164*u.u): + return (m * ac.mu0 * mu**2) / (8 * np.pi * ac.hbar**2) + +def scatteringCrossSection(B): + return 8 * np.pi * scatteringLength(B)**2 + ((32*np.pi)/45) * dipolarLength()**2 + +def calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N, T, B): #For a 3D Harmonic Trap + return (particleDensity(w_x, w_z, Power, Polarizability, N, T) * scatteringCrossSection(B) * meanThermalVelocity(T) / (2 * np.sqrt(2))).decompose() + +def calculatePSD(w_x, w_z, Power, Polarizability, N, T): + return (particleDensity(w_x, w_z, Power, Polarizability, N, T, m = 164*u.u) * thermaldeBroglieWavelength(T)**3).decompose() + +def computeTrapPotential(w_x, w_z, Power, Polarizability, options): + + axis = options['axis'] + extent = options['extent'] + gravity = options['gravity'] + astigmatism = options['astigmatism'] + + TrappingPotential = [] + TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability) + IdealTrapDepthInKelvin = (TrapDepth/ac.k_B).to(u.uK) + + projection_axis = np.array([0, 1, 0]) # default + + if axis == 0: + projection_axis = np.array([1, 0, 0]) # radial direction (X-axis) + + elif axis == 1: + projection_axis = np.array([0, 1, 0]) # propagation direction (Y-axis) + + elif axis == 2: + projection_axis = np.array([0, 0, 1]) # vertical direction (Z-axis) + + x_Positions = np.arange(-extent, extent, 1)*u.um + y_Positions = np.arange(-extent, extent, 1)*u.um + z_Positions = np.arange(-extent, extent, 1)*u.um + Positions = np.vstack((x_Positions, y_Positions, z_Positions)) * projection_axis[:, np.newaxis] + + IdealTrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability) + IdealTrappingPotential = IdealTrappingPotential * (np.ones((3, len(IdealTrappingPotential))) * projection_axis[:, np.newaxis]) + IdealTrappingPotential = (IdealTrappingPotential/ac.k_B).to(u.uK) + + if gravity and not astigmatism: + # Influence of Gravity + m = 164*u.u + gravity_axis = np.array([0, 0, 1]) + tilt_gravity = options['tilt_gravity'] + theta = options['theta'] + tilt_axis = options['tilt_axis'] + if tilt_gravity: + R = rotation_matrix(tilt_axis, np.radians(theta)) + gravity_axis = np.dot(R, gravity_axis) + gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis] + TrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability) + 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) + + elif not gravity and astigmatism: + # Influence of Astigmatism + disp_foci = options['disp_foci'] + TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci, alpha = Polarizability) + TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) + + elif gravity and astigmatism: + # Influence of Gravity and Astigmatism + m = 164*u.u + gravity_axis = np.array([0, 0, 1]) + tilt_gravity = options['tilt_gravity'] + theta = options['theta'] + tilt_axis = options['tilt_axis'] + disp_foci = options['disp_foci'] + if tilt_gravity: + R = rotation_matrix(tilt_axis, np.radians(theta)) + gravity_axis = np.dot(R, gravity_axis) + 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, alpha = Polarizability) + 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) + + 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]) + + TrapDepthsInKelvin = [IdealTrapDepthInKelvin, EffectiveTrapDepthInKelvin] + + 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') + CalculatedTrapFrequencies = [v_x, v_y, v_z] + + v, dv, popt, pcov = extractTrapFrequency(Positions, IdealTrappingPotential, IdealTrapDepthInKelvin, axis) + IdealTrapFrequency = [v, dv] + v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, EffectiveTrapDepthInKelvin, axis) + TrapFrequency = [v, dv] + ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency] + + # Gamma_elastic = calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N = 1.13 * 1e7, T = 22 * u.uK, B = 0 * u.G) + # PSD = calculatePSD(w_x, w_z, Power, Polarizability, N = 1.13 * 1e7, T = 22 * u.uK).decompose() + + return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies + def plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, pcov): v = popt[0] dv = pcov[0][0]**0.5 @@ -99,31 +230,51 @@ def plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, plt.legend(prop={'size': 12, 'weight': 'bold'}) plt.show() -def plotPotential(Positions, Powers, ComputedPotentials, axis, TrapDepthLabels): - - ## plot of the measured parameter vs. scan parameter - plt.figure(figsize=(9, 7)) - j = 0 - for i in range(np.size(ComputedPotentials, 0)): - v, dv, popt, pcov = extractTrapFrequency(Positions, ComputedPotentials[i], TrapDepthInKelvin, axis) +def generate_label(v, dv): + unit = 'Hz' + if v <= 0.0: + v = np.nan + dv = np.nan unit = 'Hz' - if v <= 0.0: - v = np.nan - dv = np.nan - unit = 'Hz' - elif v > 0.0 and orderOfMagnitude(v) > 2: - v = v / 1e3 # in kHz - dv = dv / 1e3 # in kHz - unit = 'kHz' - tf_label = '\u03BD = %.1f \u00B1 %.2f %s'% tuple([v,dv,unit]) - if np.size(ComputedPotentials, 0) == len(Powers): - plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'P = ' + str(Powers[i]) + ' W; ' + TrapDepthLabels[i] + '; ' + tf_label) + elif v > 0.0 and orderOfMagnitude(v) > 2: + v = v / 1e3 # in kHz + dv = dv / 1e3 # in kHz + unit = 'kHz' + tf_label = '\u03BD = %.1f \u00B1 %.2f %s'% tuple([v,dv,unit]) + return tf_label + +def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterateOver = [], save = False): + + plt.figure(figsize=(9, 7)) + for i in range(np.size(ComputedPotentials, 0)): + + if i % 2 == 0: + j = int(i / 2) else: - if i % 2 == 0 and j < len(Powers): - plt.plot(Positions[axis], ComputedPotentials[i][axis], '--',label = 'P = ' + str(Powers[j]) + ' W; ' + TrapDepthLabels[j] + '; ' + tf_label) - elif i % 2 != 0 and j < len(Powers): - plt.plot(Positions[axis], ComputedPotentials[i][axis], label = 'P = ' + str(Powers[j]) + ' W; ' + tf_label) - j = j + 1 + j = int((i - 1) / 2) + + IdealTrapDepthInKelvin = Params[j][0][0] + EffectiveTrapDepthInKelvin = Params[j][0][1] + + 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 listToIterateOver: + if np.size(ComputedPotentials, 0) == len(listToIterateOver): + plt.plot(Positions[axis], ComputedPotentials[i][axis], label = '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)) + 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)) + 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)) + 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)) if axis == 0: dir = 'X' elif axis == 1: @@ -135,104 +286,47 @@ def plotPotential(Positions, Powers, ComputedPotentials, axis, TrapDepthLabels): plt.tight_layout() plt.grid(visible=1) plt.legend(loc=3, prop={'size': 12, 'weight': 'bold'}) + if save: + plt.savefig('pot_' + dir + '.png') plt.show() - # plt.savefig('pot_' + dir + '.png') -if __name__ == '__main__': +if __name__ == '__main__': - # Powers = [0.1, 0.5, 2] - # Powers = [5, 20, 40] - Powers = [40] + Power = 40*u.W Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability - w_x, w_z = 34*u.um, 27.5*u.um # Beam Waists in the x and y directions - # w_x, w_z = 30*u.um, 30*u.um # Beam Waists in the x and y directions - # w_x, w_z = 20.5*u.um, 20.5*u.um - - 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 - - TrappingPotential = [] - ComputedPotentials = [] - TrapDepthLabels = [] - - gravity = True - astigmatism = True - - tilt_gravity = True - theta = 5 # in degrees - tilt_axis = [1, 0, 0] # lab space coordinates are rotated about x-axis in reference frame of beam - - 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) - - for p in Powers: - - Power = p*u.W # Single Beam Power - - TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability) - TrapDepthInKelvin = (TrapDepth/ac.k_B).to(u.uK) - TrapDepthLabels.append("Trap Depth = " + str(round(TrapDepthInKelvin.value, 2)) + " " + str(TrapDepthInKelvin.unit)) - - projection_axis = np.array([0, 1, 0]) # default - - if axis == 0: - projection_axis = np.array([1, 0, 0]) # radial direction (X-axis) - - elif axis == 1: - projection_axis = np.array([0, 1, 0]) # propagation direction (Y-axis) - - elif axis == 2: - projection_axis = np.array([0, 0, 1]) # vertical direction (Z-axis) - - x_Positions = np.arange(-extent, extent, 1)*u.um - y_Positions = np.arange(-extent, extent, 1)*u.um - z_Positions = np.arange(-extent, extent, 1)*u.um - Positions = np.vstack((x_Positions, y_Positions, z_Positions)) * projection_axis[:, np.newaxis] - - IdealTrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability) - IdealTrappingPotential = IdealTrappingPotential + np.zeros((3, len(IdealTrappingPotential))) * IdealTrappingPotential.unit - IdealTrappingPotential = (IdealTrappingPotential/ac.k_B).to(u.uK) - - if gravity and not astigmatism: - ComputedPotentials.append(IdealTrappingPotential) - # Influence of Gravity - m = 164*u.u - gravity_axis = np.array([0, 0, 1]) - if tilt_gravity: - R = rotation_matrix(tilt_axis, np.radians(theta)) - gravity_axis = np.dot(R, gravity_axis) - gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis] - TrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability) + gravitational_potential(gravity_axis_positions, m) - TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) - - elif not gravity and astigmatism: - ComputedPotentials.append(IdealTrappingPotential) - # Influence of Astigmatism - TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci, alpha = Polarizability) - TrappingPotential = TrappingPotential + np.zeros((3, len(TrappingPotential))) * TrappingPotential.unit - TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) - - elif gravity and astigmatism: - ComputedPotentials.append(IdealTrappingPotential) - # Influence of Gravity and Astigmatism - m = 164*u.u - gravity_axis = np.array([0, 0, 1]) - if tilt_gravity: - R = rotation_matrix(tilt_axis, np.radians(theta)) - gravity_axis = np.dot(R, gravity_axis) - 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, alpha = Polarizability) + gravitational_potential(gravity_axis_positions, m) - TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK) - - else: - TrappingPotential = IdealTrappingPotential - - # v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x') - # v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'z') - - # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, axis) - # plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, pcov) - - ComputedPotentials.append(TrappingPotential) + 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 + 'gravity': True, + 'astigmatism': False, + 'tilt_gravity': True, + 'theta': 5, # in degrees + 'tilt_axis': [1, 0, 0], # lab space coordinates are rotated about x-axis in reference frame of beam + '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) + } + + 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, Powers, ComputedPotentials, axis, TrapDepthLabels) \ No newline at end of file + plotPotential(Positions, ComputedPotentials, options['axis'], Params) + + # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, options['axis']) + # plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, options['axis'], popt, pcov) + + # 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['axis'], Params) \ No newline at end of file