diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py index f7fc8b8..b829611 100644 --- a/calculateDipoleTrapPotential.py +++ b/calculateDipoleTrapPotential.py @@ -4,6 +4,10 @@ import matplotlib.pyplot as plt from scipy.optimize import curve_fit from astropy import units as u, constants as ac +##################################################################### +# HELPER FUNCTIONS # +##################################################################### + def orderOfMagnitude(number): return math.floor(math.log(number, 10)) @@ -27,6 +31,10 @@ def rotation_matrix(axis, theta): [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)], [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) +##################################################################### +# BEAM PARAMETERS # +##################################################################### + # Rayleigh range def z_R(w_0:np.ndarray, lamb:float)->np.ndarray: return np.pi*w_0**2/lamb @@ -35,57 +43,9 @@ def z_R(w_0:np.ndarray, lamb:float)->np.ndarray: def w(pos, w_0, lamb): return w_0*np.sqrt(1+(pos / z_R(w_0, lamb))**2) -def trap_depth(w_1:"float|u.quantity.Quantity", w_2:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity", alpha:float)->"float|u.quantity.Quantity": - 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 gravitational_potential(positions: "np.ndarray|u.quantity.Quantity", m:"float|u.quantity.Quantity"): - return m * ac.g0 * positions - -def single_gaussian_beam_potential(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", alpha:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: - 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 = - 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 - -def astigmatic_single_gaussian_beam_potential(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", del_y:"float|u.quantity.Quantity", alpha:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: - 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 = - 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 - -def single_gaussian_beam_potential_harmonic_approximation(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", depth:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: - U = - depth * (1 - (2 * (positions[0,:]/waists[0])**2) - (2 * (positions[2,:]/waists[1])**2) - (0.5 * positions[1,:]**2 * np.sum(np.reciprocal(z_R(waists, wavelength)))**2)) - return U - -def harmonic_potential(pos, v, offset, m = 164*u.u): - U_Harmonic = ((0.5 * m * (2 * np.pi * v*u.Hz)**2 * (pos*u.um)**2)/ac.k_B).to(u.uK) + offset*u.uK - return U_Harmonic.value - -def calculateTrapFrequency(w_x, w_z, Power, Polarizability, m = 164*u.u, dir = 'x'): - TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability) - 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 - -def extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, axis): - tmp_pos = Positions[axis, :] - center_idx = np.where(tmp_pos == 0)[0][0] - lb = int(round(center_idx - len(tmp_pos)/20, 1)) - ub = int(round(center_idx + len(tmp_pos)/20, 1)) - xdata = tmp_pos[lb:ub] - tmp_pot = TrappingPotential[axis] - Potential = tmp_pot[lb:ub] - p0=[1e3, -TrapDepthInKelvin.value] - popt, pcov = curve_fit(harmonic_potential, xdata, Potential, p0) - v = popt[0] - dv = pcov[0][0]**0.5 - return v, dv, popt, pcov +##################################################################### +# RELEVANT PARAMETERS FOR EVAPORATIVE COOLING # +##################################################################### def meanThermalVelocity(T, m = 164*u.u): return 4 * np.sqrt((ac.k_B * T) /(np.pi * m)) @@ -120,6 +80,66 @@ def calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N, T, B): #Fo 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() +##################################################################### +# POTENTIALS # +##################################################################### + +def gravitational_potential(positions: "np.ndarray|u.quantity.Quantity", m:"float|u.quantity.Quantity"): + return m * ac.g0 * positions + +def single_gaussian_beam_potential(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", alpha:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: + 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 = - 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 + +def astigmatic_single_gaussian_beam_potential(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", del_y:"float|u.quantity.Quantity", alpha:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: + 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 = - 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 + +def single_gaussian_beam_potential_harmonic_approximation(positions: "np.ndarray|u.quantity.Quantity", waists: "np.ndarray|u.quantity.Quantity", depth:"float|u.quantity.Quantity"=1, wavelength:"float|u.quantity.Quantity"=1.064*u.um)->np.ndarray: + U = - depth * (1 - (2 * (positions[0,:]/waists[0])**2) - (2 * (positions[2,:]/waists[1])**2) - (0.5 * positions[1,:]**2 * np.sum(np.reciprocal(z_R(waists, wavelength)))**2)) + return U + +def harmonic_potential(pos, v, xoffset, yoffset, m = 164*u.u): + 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 + +##################################################################### +# COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS # +##################################################################### + +def trap_depth(w_1:"float|u.quantity.Quantity", w_2:"float|u.quantity.Quantity", P:"float|u.quantity.Quantity", alpha:float)->"float|u.quantity.Quantity": + 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'): + TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability) + 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 + +def extractTrapFrequency(Positions, TrappingPotential, axis): + tmp_pos = Positions[axis, :] + tmp_pot = TrappingPotential[axis] + center_idx = np.argmin(tmp_pot) + lb = int(round(center_idx - len(tmp_pot)/150, 1)) + ub = int(round(center_idx + len(tmp_pot)/150, 1)) + xdata = tmp_pos[lb:ub] + Potential = tmp_pot[lb:ub] + p0=[1e3, tmp_pos[center_idx].value, np.argmin(tmp_pot.value)] + popt, pcov = curve_fit(harmonic_potential, xdata, Potential, p0) + v = popt[0] + dv = pcov[0][0]**0.5 + return v, dv, popt, pcov + def computeTrapPotential(w_x, w_z, Power, Polarizability, options): axis = options['axis'] @@ -204,9 +224,9 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): 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) + v, dv, popt, pcov = extractTrapFrequency(Positions, IdealTrappingPotential, axis) IdealTrapFrequency = [v, dv] - v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, EffectiveTrapDepthInKelvin, axis) + v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, axis) TrapFrequency = [v, dv] ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency] @@ -215,7 +235,11 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies -def plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, pcov): +##################################################################### +# PLOT TRAP POTENTIALS # +##################################################################### + +def plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, axis, popt, pcov): v = popt[0] dv = pcov[0][0]**0.5 happrox = harmonic_potential(Positions[axis, :].value, *popt) @@ -224,7 +248,7 @@ def plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, axis, popt, plt.plot(Positions[axis, :], TrappingPotential[axis], 'ob', label = 'Gaussian Potential') plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold') plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold') - plt.ylim([-TrapDepthInKelvin.value, max(TrappingPotential[axis].value)]) + plt.ylim([-TrapDepthsInKelvin[0].value, max(TrappingPotential[axis].value)]) plt.tight_layout() plt.grid(visible=1) plt.legend(prop={'size': 12, 'weight': 'bold'}) @@ -290,6 +314,10 @@ def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterat plt.savefig('pot_' + dir + '.png') plt.show() +##################################################################### +# FUNCTION CALLS BELOW # +##################################################################### + if __name__ == '__main__': Power = 40*u.W @@ -318,8 +346,8 @@ if __name__ == '__main__': ComputedPotentials = np.asarray(ComputedPotentials) plotPotential(Positions, ComputedPotentials, options['axis'], Params) - # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, TrapDepthInKelvin, options['axis']) - # plotHarmonicFit(Positions, TrappingPotential, TrapDepthInKelvin, options['axis'], popt, pcov) + # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, options['axis']) + # plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, options['axis'], popt, pcov) # Power = [10, 20, 25, 30, 35, 40]*u.W # Single Beam Power # for p in Power: