From ac5aeb4a3ee0586d7f42be54d79efdf13f0fc23f Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 22 Feb 2023 12:42:54 +0100 Subject: [PATCH] MAJOR REWRITE --- calculateDipoleTrapPotential.py | 774 ++++++++++++++++++++++++-------- 1 file changed, 595 insertions(+), 179 deletions(-) diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py index b88510f..344753c 100644 --- a/calculateDipoleTrapPotential.py +++ b/calculateDipoleTrapPotential.py @@ -1,7 +1,8 @@ import math import numpy as np import matplotlib.pyplot as plt -from scipy import signal +import matplotlib.ticker as mtick +from scipy import signal, interpolate from scipy.optimize import curve_fit from astropy import units as u, constants as ac @@ -43,6 +44,8 @@ def modulation_function(mod_amp, n_points, func = 'arccos'): phi = np.linspace(0, 2*np.pi, n_points) mod_func = mod_amp * np.sin(phi) elif func == 'arccos': + # phi = np.linspace(0, 2*np.pi, n_points) + # mod_func = mod_amp * (2/np.pi * np.arccos(phi/np.pi-1) - 1) phi = np.linspace(0, 2*np.pi, int(n_points/2)) tmp_1 = 2/np.pi * np.arccos(phi/np.pi-1) - 1 tmp_2 = np.flip(tmp_1) @@ -80,12 +83,62 @@ def w(pos, w_0, lamb): 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 particleDensity(w_x, w_z, Power, Polarizability, N, T, m = 164*u.u, use_measured_tf = False): # For a thermal cloud + if not use_measured_tf: + 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) + 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))) + + 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))] + 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))] + + 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))) + + 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_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))) + sorted_fin_mod_dep, sorted_avg_dT = zip(*sorted(zip(fin_mod_dep, avg_dT))) + + pd = np.zeros(len(fin_mod_dep)) + dpd = np.zeros(len(fin_mod_dep)) + + for i in range(len(fin_mod_dep)): + 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() + 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 + + pd = pd*particle_density.unit + dpd = dpd*particle_density.unit + + return pd, dpd, sorted_avg_T, sorted_avg_dT, new_aspect_ratio, sorted_fin_mod_dep + def thermaldeBroglieWavelength(T, m = 164*u.u): return np.sqrt((2*np.pi*ac.hbar**2)/(m*ac.k_B*T)) @@ -141,6 +194,39 @@ 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() +def convert_modulation_depth_to_alpha(modulation_depth): + 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] + fy = [2.746, 1.278, 1.719, 1.058, 0.923, 0.994, 0.911, 1.157, 1.446, 2.191, 2.643] + dfy = [0.014, 0.007, 0.009, 0.007, 0.005, 0.004, 0.004, 0.005, 0.007, 0.009, 0.033] + + 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))] + alpha_y = [(fy[0]/y)**2 for y in fy] + dalpha_y = [alpha_y[i] * np.sqrt((dfy[0]/fy[0])**2 + (dfy[i]/fy[i])**2) for i in range(len(fy))] + + avg_alpha = [(g + h) / 2 for g, h in zip(alpha_x, alpha_y)] + sorted_fin_mod_dep, sorted_avg_alpha = zip(*sorted(zip(fin_mod_dep, avg_alpha))) + + f = interpolate.interp1d(sorted_fin_mod_dep, sorted_avg_alpha) + + return f(modulation_depth), fin_mod_dep, alpha_x, alpha_y, dalpha_x, dalpha_y + +def convert_modulation_depth_to_temperature(modulation_depth): + fin_mod_dep = [1.0, 0.8, 0.6, 0.4, 0.2, 0.0, 0.9, 0.7, 0.5, 0.3, 0.1] + T_x = [22.1, 27.9, 31.7, 42.2, 98.8, 145.8, 27.9, 33.8, 42.4, 61.9, 136.1] + dT_x = [1.7, 2.6, 2.4, 3.7, 1.1, 0.6, 2.2, 3.2, 1.7, 2.2, 1.2] + T_y = [13.13, 14.75, 18.44, 26.31, 52.55, 92.9, 13.54, 16.11, 21.15, 35.81, 85.8] + dT_y = [0.05, 0.05, 0.07, 0.16, 0.28, 0.7, 0.04, 0.07, 0.10, 0.21, 0.8] + + avg_T = [(g + h) / 2 for g, h in zip(T_x, T_y)] + sorted_fin_mod_dep, sorted_avg_T = zip(*sorted(zip(fin_mod_dep, avg_T))) + + f = interpolate.interp1d(sorted_fin_mod_dep, sorted_avg_T) + + return f(modulation_depth), fin_mod_dep, T_x, T_y, dT_x, dT_y + ##################################################################### # POTENTIALS # ##################################################################### @@ -180,6 +266,23 @@ def gaussian_potential(pos, amp, waist, xoffset, yoffset): U_Gaussian = amp * np.exp(-2 * ((pos + xoffset) / waist)**2) + yoffset return U_Gaussian +def crossed_beam_potential(positions, theta, waists, P, alpha, wavelength=1.064*u.um): + + beam_1_positions = positions + A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:], waists[0][0], wavelength)*w(beam_1_positions[1,:], waists[0][1], wavelength)) + U_1_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) + U_1 = - U_1_tilde * A_1 * np.exp(-2 * ((beam_1_positions[0,:]/w(beam_1_positions[1,:], waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:], waists[0][1], wavelength))**2)) + + R = rotation_matrix([0, 0, 1], np.radians(theta)) + beam_2_positions = np.dot(R, beam_1_positions) + A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:], waists[1][0], wavelength)*w(beam_2_positions[1,:], waists[1][1], wavelength)) + U_2_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) + U_2 = - U_2_tilde * A_2 * np.exp(-2 * ((beam_2_positions[0,:]/w(beam_2_positions[1,:], waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:], waists[1][1], wavelength))**2)) + + U = U_1 + U_2 + + return U + ##################################################################### # COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS # ##################################################################### @@ -220,6 +323,7 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): gravity = options['gravity'] astigmatism = options['astigmatism'] modulation = options['modulation'] + crossed = options['crossed'] if modulation: aspect_ratio = options['aspect_ratio'] @@ -241,14 +345,25 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): elif axis == 2: projection_axis = np.array([0, 0, 1]) # vertical direction (Z-axis) + else: + projection_axis = np.array([1, 1, 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 not crossed: + 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) + + else: + theta = options['theta'] + 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 = 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 @@ -291,27 +406,31 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): 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]) + if not crossed: + 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]) + else: + EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin + + 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, axis) + IdealTrapFrequency = [v, dv] + v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, axis) + TrapFrequency = [v, dv] + ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency] + + return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies + else: - EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin - - 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, axis) - IdealTrapFrequency = [v, dv] - v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, axis) - TrapFrequency = [v, dv] - ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency] - - return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies + return TrappingPotential def extractWaist(Positions, TrappingPotential): tmp_pos = Positions.value @@ -321,8 +440,8 @@ def extractWaist(Positions, TrappingPotential): TrapMinimum = tmp_pot[center_idx] TrapCenter = tmp_pos[center_idx] - lb = int(round(center_idx - len(tmp_pot)/10, 1)) - ub = int(round(center_idx + len(tmp_pot)/10, 1)) + lb = int(round(center_idx - len(tmp_pot)/30, 1)) + ub = int(round(center_idx + len(tmp_pot)/30, 1)) xdata = tmp_pos[lb:ub] Potential = tmp_pot[lb:ub] @@ -330,10 +449,87 @@ def extractWaist(Positions, TrappingPotential): popt, pcov = curve_fit(gaussian_potential, xdata, Potential, p0) return popt, pcov +def computeIntensityProfileAndPotentials(Power, waists, alpha, wavelength, options): + w_x = waists[0] + w_z = waists[1] + extent = options['extent'] + modulation = options['modulation'] + mod_func = options['modulation_function'] + + if not modulation: + extent = 50 + 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 + + idx = np.where(y_Positions==0)[0][0] + + alpha = Polarizability + wavelength = 1.064*u.um + + xm,ym,zm = np.meshgrid(x_Positions, y_Positions, z_Positions, sparse=True, indexing='ij') + + ## Single Gaussian Beam + A = 2*Power/(np.pi*w(ym, w_x, wavelength)*w(ym, w_z, wavelength)) + intensity_profile = A * np.exp(-2 * ((xm/w(ym, w_x, wavelength))**2 + (zm/w(ym, w_z, wavelength))**2)) + I = intensity_profile[:, idx, :].to(u.MW/(u.cm*u.cm)) + + U_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) + U = - U_tilde * I + U = (U/ac.k_B).to(u.uK) + + return [x_Positions, z_Positions], [w_x.value, 0, w_z.value, 0], I, U, [0, 0, 0, 0] + + else: + mod_amp = options['modulation_amplitude'] + 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 + + mod_amp = mod_amp * w_x + n_points = len(x_Positions) + dx, xmod_Positions = modulation_function(mod_amp, n_points, func = mod_func) + + idx = np.where(y_Positions==0)[0][0] + + xm,ym,zm,xmodm = np.meshgrid(x_Positions, y_Positions, z_Positions, xmod_Positions, sparse=True, indexing='ij') + + ## Single Modulated Gaussian Beam + A = 2*Power/(np.pi*w(y_Positions[idx] , w_x, wavelength)*w(y_Positions[idx], w_z, wavelength)) + intensity_profile = A * 1/(2*mod_amp) * np.trapz(np.exp(-2 * (((xmodm - xm)/w(ym, w_x, wavelength))**2 + (zm/w(ym, w_z, wavelength))**2)), dx = dx, axis = -1) + I = intensity_profile[:, idx, :].to(u.MW/(u.cm*u.cm)) + + U_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) + U = - U_tilde * I + U = (U/ac.k_B).to(u.uK) + + poptx, pcovx = extractWaist(x_Positions, U[:, np.where(z_Positions==0)[0][0]]) + poptz, pcovz = extractWaist(z_Positions, U[np.where(x_Positions==0)[0][0], :]) + + extracted_waist_x = poptx[1] + dextracted_waist_x = pcovx[1][1]**0.5 + extracted_waist_z = poptz[1] + dextracted_waist_z = pcovz[1][1]**0.5 + + return [x_Positions, z_Positions], [extracted_waist_x, dextracted_waist_x, extracted_waist_z, dextracted_waist_z], I, U, [poptx, pcovx, poptz, pcovz] + ##################################################################### # PLOTTING # ##################################################################### +def generate_label(v, dv): + 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]) + return tf_label + def plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, axis, popt, pcov): v = popt[0] dv = pcov[0][0]**0.5 @@ -379,26 +575,13 @@ def plotGaussianFit(Positions, TrappingPotential, popt, pcov): bx.set_title('Fit Residuals') plt.plot(Positions, TrappingPotential - gapprox, 'ob') plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold') - plt.ylabel('$U_{trap} - U_{Harmonic}$', fontsize= 12, fontweight='bold') + plt.ylabel('$U_{trap} - U_{Gaussian}$', fontsize= 12, fontweight='bold') plt.xlim([-10, 10]) plt.ylim([-1, 1]) plt.grid(visible=1) plt.tight_layout() plt.show() -def generate_label(v, dv): - 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]) - return tf_label - def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterateOver = [], save = False): plt.figure(figsize=(9, 7)) @@ -448,108 +631,195 @@ def plotPotential(Positions, ComputedPotentials, axis, Params = [], listToIterat plt.savefig('pot_' + dir + '.png') plt.show() -def plotIntensityProfileAndPotentials(Power, waists, alpha, wavelength, options): - w_x = waists[0] - w_z = waists[1] - extent = options['extent'] - modulation = options['modulation'] - mod_func = options['modulation_function'] - - if not modulation: - extent = 50 - 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 - - idx = np.where(y_Positions==0)[0][0] +def plotIntensityProfileAndPotentials(positions, waists, I, U): + + x_Positions = positions[0] + z_Positions = positions[1] + + w_x = waists[0] + dw_x = waists[1] + w_z = waists[2] + dw_x = waists[3] - alpha = Polarizability - wavelength = 1.064*u.um - - xm,ym,zm = np.meshgrid(x_Positions, y_Positions, z_Positions, sparse=True, indexing='ij') + ar = w_x/w_z + dar = ar * np.sqrt((dw_x/w_x)**2 + (dw_x/w_z)**2) + + fig = plt.figure(figsize=(12, 6)) + ax = fig.add_subplot(121) + ax.set_title('Intensity Profile ($MW/cm^2$)\n Aspect Ratio = %.2f \u00B1 %.2f um'% tuple([ar,dar])) + im = plt.imshow(np.transpose(I.value), cmap="coolwarm", extent=[np.min(x_Positions.value), np.max(x_Positions.value), np.min(z_Positions.value), np.max(z_Positions.value)]) + plt.xlabel('X - Horizontal (um)', fontsize= 12, fontweight='bold') + plt.ylabel('Z - Vertical (um)', fontsize= 12, fontweight='bold') + ax.set_aspect('equal') + fig.colorbar(im, fraction=0.046, pad=0.04, orientation='vertical') + + bx = fig.add_subplot(122) + bx.set_title('Trap Potential') + plt.plot(x_Positions, U[:, np.where(z_Positions==0)[0][0]], label = 'X - Horizontal') + plt.plot(z_Positions, U[np.where(x_Positions==0)[0][0], :], label = 'Z - Vertical') + plt.ylim(top = 0) + plt.xlabel('Extent (um)', fontsize= 12, fontweight='bold') + plt.ylabel('Depth (uK)', fontsize= 12, fontweight='bold') + plt.tight_layout() + plt.grid(visible=1) + plt.legend(prop={'size': 12, 'weight': 'bold'}) + plt.show() - ## Single Gaussian Beam - A = 2*Power/(np.pi*w(ym, w_x, wavelength)*w(ym, w_z, wavelength)) - I = A * np.exp(-2 * ((xm/w(ym, w_x, wavelength))**2 + (zm/w(ym, w_z, wavelength))**2)) - I = np.transpose(I.to(u.MW/(u.cm*u.cm))) - - U_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) - U = - U_tilde * I - U = (U/ac.k_B).to(u.uK) +def plotAlphas(): - fig = plt.figure(figsize=(12, 6)) - ax = fig.add_subplot(121) - ax.set_title('Intensity Profile ($MW/cm^2$)\n Aspect Ratio = %.2f' %(w_x/w_z)) - im = plt.imshow(I[:,idx,:].value, cmap="coolwarm", extent=[np.min(x_Positions.value), np.max(x_Positions.value), np.min(z_Positions.value), np.max(z_Positions.value)]) - plt.xlabel('X - Horizontal (um)', fontsize= 12, fontweight='bold') - plt.ylabel('Z - Vertical (um)', fontsize= 12, fontweight='bold') - ax.set_aspect('equal') - fig.colorbar(im, fraction=0.046, pad=0.04, orientation='vertical') - - bx = fig.add_subplot(122) - bx.set_title('Trap Potential') - plt.plot(x_Positions, U[np.where(x_Positions==0)[0][0], idx, :], label = 'X - Horizontal') - plt.plot(z_Positions, U[:, idx, np.where(z_Positions==0)[0][0]], label = 'Z - Vertical') - plt.ylim(top = 0) - plt.xlabel('Extent (um)', fontsize= 12, fontweight='bold') - plt.ylabel('Depth (uK)', fontsize= 12, fontweight='bold') - plt.tight_layout() - plt.grid(visible=1) - plt.legend(prop={'size': 12, 'weight': 'bold'}) - plt.show() - else: - mod_amp = options['modulation_amplitude'] - 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 + modulation_depth = np.arange(0, 1.1, 0.1) + Alphas, fin_mod_dep, alpha_x, alpha_y, dalpha_x, dalpha_y = convert_modulation_depth_to_alpha(modulation_depth) - mod_amp = mod_amp * w_x - n_points = len(x_Positions) - dx, xmod_Positions = modulation_function(mod_amp, n_points, func = mod_func) - - idx = np.where(y_Positions==0)[0][0] + plt.figure() + plt.errorbar(fin_mod_dep, alpha_x, yerr = dalpha_x, fmt= 'ob', markersize=5, capsize=5) + plt.errorbar(fin_mod_dep, alpha_y, yerr = dalpha_y, fmt= 'or', markersize=5, capsize=5) + plt.plot(modulation_depth, Alphas, '--g') + plt.xlabel('Modulation depth', fontsize= 12, fontweight='bold') + plt.ylabel('$\\alpha$', fontsize= 12, fontweight='bold') + plt.tight_layout() + plt.grid(visible=1) + plt.show() - xm,ym,zm,xmodm = np.meshgrid(x_Positions, y_Positions, z_Positions, xmod_Positions, sparse=True, indexing='ij') +def plotTemperatures(w_x, w_z, plot_against_mod_depth = True): - ## Single Modulated Gaussian Beam - A = 2*Power/(np.pi*w(y_Positions[idx] , w_x, wavelength)*w(y_Positions[idx], w_z, wavelength)) - intensity_profile = A * 1/(2*mod_amp) * np.trapz(np.exp(-2 * (((xmodm - xm)/w(ym, w_x, wavelength))**2 + (zm/w(ym, w_z, wavelength))**2)), dx = dx, axis = -1) - I = intensity_profile[:, idx, :].to(u.MW/(u.cm*u.cm)) - - U_tilde = (1 / (2 * ac.eps0 * ac.c)) * alpha * (4 * np.pi * ac.eps0 * ac.a0**3) - U = - U_tilde * I - U = (U/ac.k_B).to(u.uK) + modulation_depth = np.arange(0, 1.1, 0.1) + w_xs = w_x * convert_modulation_depth_to_alpha(modulation_depth)[0] + new_aspect_ratio = w_xs / w_z + Temperatures, fin_mod_dep, T_x, T_y, dT_x, dT_y = convert_modulation_depth_to_temperature(modulation_depth) + measured_aspect_ratio = (w_x * convert_modulation_depth_to_alpha(fin_mod_dep)[0]) / w_z - poptx, pcovx = extractWaist(x_Positions, U[:, np.where(z_Positions==0)[0][0]]) - poptz, pcovz = extractWaist(z_Positions, U[np.where(x_Positions==0)[0][0], :]) - - extracted_waist_x = poptx[1] - dextracted_waist_x = pcovx[1][1]**0.5 - extracted_waist_z = poptz[1] - dextracted_waist_z = pcovz[1][1]**0.5 - ar = extracted_waist_x/extracted_waist_z - dar = ar * np.sqrt((dextracted_waist_x/extracted_waist_x)**2 + (dextracted_waist_z/extracted_waist_z)**2) + plt.figure() + if plot_against_mod_depth: + plt.errorbar(fin_mod_dep, T_x, yerr = dT_x, fmt= 'ob', markersize=5, capsize=5) + plt.errorbar(fin_mod_dep, T_y, yerr = dT_y, fmt= 'or', markersize=5, capsize=5) + plt.plot(modulation_depth, Temperatures, '--g') + xlabel = 'Modulation depth' + else: + plt.errorbar(measured_aspect_ratio, T_x, yerr = dT_x, fmt= 'ob', markersize=5, capsize=5) + plt.errorbar(measured_aspect_ratio, T_y, yerr = dT_y, fmt= 'or', markersize=5, capsize=5) + plt.plot(new_aspect_ratio, Temperatures, '--g') + xlabel = 'Aspect Ratio' - fig = plt.figure(figsize=(12, 6)) - ax = fig.add_subplot(121) - ax.set_title('Intensity Profile ($MW/cm^2$)\n Aspect Ratio = %.2f \u00B1 %.2f um'% tuple([ar,dar])) - im = plt.imshow(np.transpose(I.value), cmap="coolwarm", extent=[np.min(x_Positions.value), np.max(x_Positions.value), np.min(z_Positions.value), np.max(z_Positions.value)]) - plt.xlabel('X - Horizontal (um)', fontsize= 12, fontweight='bold') - plt.ylabel('Z - Vertical (um)', fontsize= 12, fontweight='bold') - ax.set_aspect('equal') - fig.colorbar(im, fraction=0.046, pad=0.04, orientation='vertical') - - bx = fig.add_subplot(122) - bx.set_title('Trap Potential') - plt.plot(x_Positions, U[:, np.where(z_Positions==0)[0][0]], label = 'X - Horizontal') - plt.plot(z_Positions, U[np.where(x_Positions==0)[0][0], :], label = 'Z - Vertical') - plt.ylim(top = 0) - plt.xlabel('Extent (um)', fontsize= 12, fontweight='bold') - plt.ylabel('Depth (uK)', fontsize= 12, fontweight='bold') - plt.tight_layout() - plt.grid(visible=1) - plt.legend(prop={'size': 12, 'weight': 'bold'}) - plt.show() + plt.xlabel(xlabel, fontsize= 12, fontweight='bold') + plt.ylabel('Temperature (uK)', fontsize= 12, fontweight='bold') + plt.tight_layout() + plt.grid(visible=1) + plt.show() + +def plotTrapFrequencies(v_x, v_y, v_z, modulation_depth, new_aspect_ratio, plot_against_mod_depth = True): + fig, ax3 = plt.subplots(figsize=(8, 6)) + + if plot_against_mod_depth: + ln1 = ax3.plot(modulation_depth, v_x, '-ob', label = 'v_x') + ln2 = ax3.plot(modulation_depth, v_z, '-^b', label = 'v_z') + ax4 = ax3.twinx() + ln3 = ax4.plot(modulation_depth, v_y, '-*r', label = 'v_y') + xlabel = 'Modulation depth' + else: + ln1 = ax3.plot(new_aspect_ratio, v_x, '-ob', label = 'v_x') + ln2 = ax3.plot(new_aspect_ratio, v_z, '-^b', label = 'v_z') + ax4 = ax3.twinx() + ln3 = ax4.plot(new_aspect_ratio, v_y, '-*r', label = 'v_y') + xlabel = 'Aspect Ratio' + + ax3.set_xlabel(xlabel, fontsize= 12, fontweight='bold') + ax3.set_ylabel('Trap Frequency (Hz)', fontsize= 12, fontweight='bold') + ax3.tick_params(axis="y", labelcolor='b') + ax4.set_ylabel('Trap Frequency (Hz)', fontsize= 12, fontweight='bold') + ax4.tick_params(axis="y", labelcolor='r') + plt.tight_layout() + plt.grid(visible=1) + lns = ln1+ln2+ln3 + labs = [l.get_label() for l in lns] + ax3.legend(lns, labs, prop={'size': 12, 'weight': 'bold'}) + plt.show() + +def plotMeasuredTrapFrequencies(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] + 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] + dalpha_y = [alpha_y[i] * np.sqrt((dfy[0]/fy[0])**2 + (dfy[i]/fy[i])**2) for i in range(len(fy))] + + avg_alpha = [(g + h) / 2 for g, h in zip(alpha_x, alpha_y)] + new_aspect_ratio = (w_x * avg_alpha) / w_z + + + if plot_against_mod_depth: + fig, ax1 = plt.subplots(figsize=(8, 6)) + ax2 = ax1.twinx() + ax1.errorbar(fin_mod_dep, 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) + ax1.errorbar(fin_mod_dep, fz, yerr = dfz, fmt= '^b', label = 'v_z', markersize=5, capsize=5) + ax1.set_xlabel('Modulation depth', fontsize= 12, fontweight='bold') + ax1.set_ylabel('Trap Frequency (kHz)', fontsize= 12, fontweight='bold') + ax1.tick_params(axis="y", labelcolor='b') + ax2.set_ylabel('Trap Frequency (Hz)', fontsize= 12, fontweight='bold') + ax2.tick_params(axis="y", labelcolor='r') + h1, l1 = ax1.get_legend_handles_labels() + h2, l2 = ax2.get_legend_handles_labels() + ax1.legend(h1+h2, l1+l2, loc=0) + else: + plt.figure() + plt.errorbar(new_aspect_ratio, fx, yerr = dfx, fmt= 'or', label = 'v_x', markersize=5, capsize=5) + plt.errorbar(new_aspect_ratio, fz, yerr = dfz, fmt= '^b', label = 'v_z', markersize=5, capsize=5) + plt.xlabel('Aspect Ratio', fontsize= 12, fontweight='bold') + plt.ylabel('Trap Frequency (kHz)', fontsize= 12, fontweight='bold') + plt.legend(prop={'size': 12, 'weight': 'bold'}) + + plt.tight_layout() + plt.grid(visible=1) + plt.show() + +def plotRatioOfTrapFrequencies(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] + 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() + + if plot_against_mod_depth: + plt.errorbar(modulation_depth, fx/v_x, yerr = dfx/v_x, fmt= 'or', label = 'b/w horz TF', markersize=5, capsize=5) + plt.errorbar(modulation_depth, fy/v_y, yerr = dfy/v_y, fmt= '*g', label = 'b/w axial TF', markersize=5, capsize=5) + plt.errorbar(modulation_depth, fz/v_z, yerr = dfz/v_z, fmt= '^b', label = 'b/w vert TF', markersize=5, capsize=5) + xlabel = 'Modulation depth' + else: + plt.errorbar(new_aspect_ratio, fx/v_x, yerr = dfx/v_x, fmt= 'or', label = 'b/w horz TF', markersize=5, capsize=5) + plt.errorbar(new_aspect_ratio, fy/v_y, yerr = dfy/v_y, fmt= '*g', label = 'b/w axial TF', markersize=5, capsize=5) + plt.errorbar(new_aspect_ratio, fz/v_z, yerr = dfz/v_z, fmt= '^b', label = 'b/w vert TF', markersize=5, capsize=5) + xlabel = 'Aspect Ratio' + + plt.xlabel(xlabel, fontsize= 12, fontweight='bold') + plt.ylabel('Ratio', fontsize= 12, fontweight='bold') + plt.tight_layout() + plt.grid(visible=1) + plt.legend(prop={'size': 12, 'weight': 'bold'}) + plt.show() def plotScatteringLengths(): BField = np.arange(0, 2.59, 1e-3) * u.G @@ -572,6 +842,30 @@ def plotScatteringLengths(): plt.grid(visible=1) plt.show() +def plotCollisionRatesAndPSD(Gamma_elastic, PSD, modulation_depth, new_aspect_ratio, plot_against_mod_depth = True): + fig, ax1 = plt.subplots(figsize=(8, 6)) + ax2 = ax1.twinx() + + if plot_against_mod_depth: + ax1.plot(modulation_depth, Gamma_elastic, '-ob') + ax2.plot(modulation_depth, PSD, '-*r') + ax2.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e')) + xlabel = 'Modulation depth' + else: + ax1.plot(new_aspect_ratio, Gamma_elastic, '-ob') + ax2.plot(new_aspect_ratio, PSD, '-*r') + ax2.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e')) + xlabel = 'Aspect Ratio' + + ax1.set_xlabel(xlabel, fontsize= 12, fontweight='bold') + ax1.set_ylabel('Elastic Collision Rate', 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() + ##################################################################### # RUN SCRIPT WITH OPTIONS BELOW # ##################################################################### @@ -587,20 +881,22 @@ if __name__ == '__main__': # 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': 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 - 'modulation': True, - 'aspect_ratio': 3.67, - '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) - } + # 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': True, + # 'aspect_ratio': 3.67, + # '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) + # } - """Plot ideal and trap potential resulting for given parameters only""" + """Plot ideal trap potential resulting for given parameters only""" # ComputedPotentials = [] # Params = [] @@ -631,50 +927,170 @@ if __name__ == '__main__': """Plot transverse intensity profile and trap potential resulting for given parameters only""" # options = { - # 'extent': 50, # range of spatial coordinates in one direction to calculate trap potential over + # 'extent': 60, # range of spatial coordinates in one direction to calculate trap potential over # 'modulation': True, # 'modulation_function': 'arccos', - # 'modulation_amplitude': 2.12 + # 'modulation_amplitude': 2.16 # } - # plotIntensityProfileAndPotentials(Power, [w_x, w_z], Polarizability, Wavelength, options) + # 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, 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 = 2.5 * u.G - modulation = True + BField = 1.4 * u.G + # modulation_depth = np.arange(0, 1.0, 0.02) + + # 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 - if modulation: - aspect_ratio = 3.67 - init_ar = w_x / w_z - w_x = w_x * (aspect_ratio / init_ar) - Temperature = 20 * u.uK - else: - Temperature = 100 * 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() + plot_against_mod_depth = True - 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)) + # # 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)) - 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') + # 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 - 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)) + """Plot alphas""" - print('a_s = %.2f ' %(scatteringLength(BField)[0] / ac.a0)) + # 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() + + """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() + + """Plot ideal crossed beam trap potential resulting for given parameters only""" + + # Powers = [40, 11] * 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 + + # options = { + # 'axis': 3, # axis referenced to the beam along which you want the dipole trap potential + # 'extent': 1e2, # range of spatial coordinates in one direction to calculate trap potential over + # 'crossed': True, + # 'theta': 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) + # } + + # 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(TrapPotential[0]) + # plt.show() + \ No newline at end of file