diff --git a/ODT Calculator/calculateDipoleTrapPotential.py b/ODT Calculator/Calculator.py similarity index 50% rename from ODT Calculator/calculateDipoleTrapPotential.py rename to ODT Calculator/Calculator.py index b5654e8..94ae4b5 100644 --- a/ODT Calculator/calculateDipoleTrapPotential.py +++ b/ODT Calculator/Calculator.py @@ -1,9 +1,6 @@ -import math import numpy as np -import matplotlib.pyplot as plt -import matplotlib.ticker as mtick -from scipy import signal, interpolate from scipy.optimize import curve_fit +from scipy import interpolate from astropy import units as u, constants as ac DY_POLARIZABILITY = 184.4 # in a.u, most precise measured value of Dy polarizability @@ -11,77 +8,7 @@ DY_MASS = 164*u.u DY_DIPOLE_MOMENT = 9.93 * ac.muB ##################################################################### -# HELPER FUNCTIONS # -##################################################################### - -def orderOfMagnitude(number): - return math.floor(math.log(number, 10)) - -def rotation_matrix(axis, theta): - """ - Return the rotation matrix associated with counterclockwise rotation about - the given axis by theta radians. - In 2-D it is just, - thetaInRadians = np.radians(theta) - c, s = np.cos(thetaInRadians), np.sin(thetaInRadians) - R = np.array(((c, -s), (s, c))) - In 3-D, one way to do it is use the Euler-Rodrigues Formula as is done here - """ - axis = np.asarray(axis) - axis = axis / math.sqrt(np.dot(axis, axis)) - a = math.cos(theta / 2.0) - b, c, d = -axis * math.sin(theta / 2.0) - aa, bb, cc, dd = a * a, b * b, c * c, d * d - bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d - return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)], - [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)], - [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) - -def find_nearest(array, value): - array = np.asarray(array) - idx = (np.abs(array - value)).argmin() - return idx - -def modulation_function(mod_amp, n_points, func = 'arccos'): - - if func == 'sin': - 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) - mod_func = mod_amp * np.concatenate((tmp_1, tmp_2)) - elif func == 'triangle': - phi = np.linspace(0, 2*np.pi, n_points) - mod_func = mod_amp * signal.sawtooth(phi, width = 0.5) # width of 0.5 gives symmetric rising triangle ramp - elif func == 'square': - phi = np.linspace(0, 1.99*np.pi, n_points) - mod_func = mod_amp * signal.square(phi, duty = 0.5) - else: - mod_func = None - - if mod_func is not None: - dx = (max(mod_func) - min(mod_func))/(2*n_points) - - return dx, mod_func - -##################################################################### -# BEAM PARAMETERS # -##################################################################### - -# Rayleigh length -def z_R(w_0, lamb)->np.ndarray: - return np.pi*w_0**2/lamb - -# Beam Radius -def w(pos, w_0, lamb): - return w_0*np.sqrt(1+(pos / z_R(w_0, lamb))**2) - -##################################################################### -# CALCULATIONS +# AUXILIARY COMPUTATIONS ##################################################################### def calculateHeatingRate(w_x, w_y, P, linewidth, detuning, wavelength = 1.064*u.um): @@ -218,104 +145,6 @@ def convert_modulation_depth_to_temperature(modulation_depth): return f(modulation_depth), fin_mod_dep, T_x, T_y, dT_x, dT_y -##################################################################### -# POTENTIALS # -##################################################################### - -def gravitational_potential(positions, m): - return m * ac.g0 * positions - -def single_gaussian_beam_potential(positions, waists, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um): - 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, waists, del_y, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um): - 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 modulated_single_gaussian_beam_potential(positions, waists, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um, mod_amp=1): - mod_amp = mod_amp * waists[0] - n_points = len(positions[0,:]) - dx, x_mod = modulation_function(mod_amp, n_points, func = 'arccos') - 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) - dU = np.zeros(2*n_points) - for i in range(len(x_mod)): - dU = np.vstack((dU, np.exp(-2 * (np.subtract(x_mod[i], positions[0,:])/w(positions[1,:], waists[0], wavelength))**2))) - U = - U_tilde * A * 1/(2*mod_amp) * np.trapz(dU, dx = dx, axis = 0) - return U - -def harmonic_potential(pos, v, xoffset, yoffset, m = DY_MASS): - 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 - -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, waists, P, options, alpha = DY_POLARIZABILITY, wavelength=1.064*u.um): - - delta = options['delta'] - - foci_shift = options['foci_shift'] - focus_shift_beam_1 = foci_shift[0] - focus_shift_beam_2 = foci_shift[1] - - beam_disp = options['beam_disp'] - beam_1_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[0])).T * beam_disp[0].unit - beam_2_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[1])).T * beam_disp[1].unit - - beam_1_positions = positions + beam_1_disp - A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][0], wavelength)*w(beam_1_positions[1,:] + focus_shift_beam_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,:] + focus_shift_beam_1, waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][1], wavelength))**2)) - - R = rotation_matrix([0, 0, 1], np.radians(delta)) - beam_2_positions = np.dot(R, positions + beam_2_disp) - A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][0], wavelength)*w(beam_2_positions[1,:] + focus_shift_beam_2, 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,:] + focus_shift_beam_2, waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][1], wavelength))**2)) - - U = U_1 + U_2 - - return U - -def astigmatic_crossed_beam_potential(positions, waists, P, options, alpha = DY_POLARIZABILITY, wavelength=1.064*u.um): - - delta = options['delta'] - - del_y = options['foci_disp_crossed'] - del_y_1 = del_y[0] - del_y_2 = del_y[1] - - - foci_shift = options['foci_shift'] - focus_shift_beam_1 = foci_shift[0] - focus_shift_beam_2 = foci_shift[1] - - beam_disp = options['beam_disp'] - beam_1_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[0])).T * beam_disp[0].unit - beam_2_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[1])).T * beam_disp[1].unit - - beam_1_positions = positions + beam_1_disp - A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] - (del_y_1/2) + focus_shift_beam_1, waists[0][0], wavelength)*w(beam_1_positions[1,:] + (del_y_1/2) + focus_shift_beam_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,:] - (del_y_1/2) + focus_shift_beam_1, waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + (del_y_1/2) + focus_shift_beam_1, waists[0][1], wavelength))**2)) - - R = rotation_matrix([0, 0, 1], np.radians(delta)) - beam_2_positions = np.dot(R, positions + beam_2_disp) - A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:] - (del_y_2/2) + focus_shift_beam_2, waists[1][0], wavelength)*w(beam_2_positions[1,:] + (del_y_2/2) + focus_shift_beam_2, 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,:] - (del_y_2/2) + focus_shift_beam_2, waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + (del_y_2/2) + focus_shift_beam_2, waists[1][1], wavelength))**2)) - - U = U_1 + U_2 - - return U - ##################################################################### # COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS # ##################################################################### @@ -676,341 +505,3 @@ def computeIntensityProfileAndPotentials(Power, waists, wavelength, options, alp 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 - happrox = harmonic_potential(Positions[axis, :].value, *popt) - fig = plt.figure(figsize=(12, 6)) - ax = fig.add_subplot(121) - ax.set_title('Fit to Potential') - plt.plot(Positions[axis, :].value, happrox, '-r', label = '\u03BD = %.1f \u00B1 %.2f Hz'% tuple([v,dv])) - 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([-TrapDepthsInKelvin[0].value, max(TrappingPotential[axis].value)]) - plt.grid(visible=1) - plt.legend(prop={'size': 12, 'weight': 'bold'}) - - bx = fig.add_subplot(122) - bx.set_title('Fit Residuals') - plt.plot(Positions[axis, :].value, TrappingPotential[axis].value - happrox, 'ob') - plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold') - plt.ylabel('$U_{trap} - U_{Harmonic}$', fontsize= 12, fontweight='bold') - plt.xlim([-10, 10]) - plt.ylim([-1e-2, 1e-2]) - plt.grid(visible=1) - plt.tight_layout() - plt.show() - -def plotGaussianFit(Positions, TrappingPotential, popt, pcov): - extracted_waist = popt[1] - dextracted_waist = pcov[1][1]**0.5 - gapprox = gaussian_potential(Positions, *popt) - fig = plt.figure(figsize=(12, 6)) - ax = fig.add_subplot(121) - ax.set_title('Fit to Potential') - plt.plot(Positions, gapprox, '-r', label = 'waist = %.1f \u00B1 %.2f um'% tuple([extracted_waist,dextracted_waist])) - plt.plot(Positions, TrappingPotential, 'ob', label = 'Gaussian Potential') - plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold') - plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold') - plt.ylim([min(TrappingPotential), max(TrappingPotential)]) - plt.grid(visible=1) - plt.legend(prop={'size': 12, 'weight': 'bold'}) - - bx = fig.add_subplot(122) - bx.set_title('Fit Residuals') - plt.plot(Positions, TrappingPotential - gapprox, 'ob') - plt.xlabel('Distance (um)', 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 plotPotential(Positions, ComputedPotentials, options, Params = [], listToIterateOver = [], save = False): - - axis = options['axis'] - - plt.figure(figsize=(9, 7)) - for i in range(np.size(ComputedPotentials, 0)): - - if i % 2 == 0: - j = int(i / 2) - else: - 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] - - 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): - 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 - Horizontal' - elif axis == 1: - dir = 'Y - Propagation' - else: - dir = 'Z - Vertical' - - plt.ylim(top = 0) - plt.xlabel(dir + ' Direction (um)', fontsize= 12, fontweight='bold') - plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold') - 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() - -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] - - 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() - -def plotAlphas(): - - 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) - - plt.figure() - plt.errorbar(fin_mod_dep, alpha_x, yerr = dalpha_x, fmt= 'ob', label = 'From Horz TF', markersize=5, capsize=5) - plt.errorbar(fin_mod_dep, alpha_y, yerr = dalpha_y, fmt= 'or', label = 'From Vert TF', 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.legend(prop={'size': 12, 'weight': 'bold'}) - plt.show() - -def plotTemperatures(w_x, w_z, plot_against_mod_depth = True): - - 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 - - plt.figure() - if plot_against_mod_depth: - plt.errorbar(fin_mod_dep, T_x, yerr = dT_x, fmt= 'ob', label = 'Horz direction', markersize=5, capsize=5) - plt.errorbar(fin_mod_dep, T_y, yerr = dT_y, fmt= 'or', label = 'Vert direction', 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', label = 'Horz direction', markersize=5, capsize=5) - plt.errorbar(measured_aspect_ratio, T_y, yerr = dT_y, fmt= 'or', label = 'Vert direction', markersize=5, capsize=5) - plt.plot(new_aspect_ratio, Temperatures, '--g') - xlabel = 'Aspect Ratio' - - plt.xlabel(xlabel, fontsize= 12, fontweight='bold') - plt.ylabel('Temperature (uK)', fontsize= 12, fontweight='bold') - plt.tight_layout() - plt.grid(visible=1) - plt.legend(prop={'size': 12, 'weight': 'bold'}) - 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, '-or', label = 'v_x') - ln2 = ax3.plot(modulation_depth, v_z, '-^b', label = 'v_z') - ax4 = ax3.twinx() - ln3 = ax4.plot(modulation_depth, v_y, '-*g', label = 'v_y') - xlabel = 'Modulation depth' - else: - ln1 = ax3.plot(new_aspect_ratio, v_x, '-or', 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, '-*g', 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='g') - 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(fx, dfx, fy, dfy, fz, dfz, modulation_depth_radial, modulation_depth_axial, w_x, w_z, plot_against_mod_depth = True): - - 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(modulation_depth_radial, fx, yerr = dfx, fmt= 'or', label = 'v_x', markersize=5, capsize=5) - ax2.errorbar(modulation_depth_axial, fy, yerr = dfy, fmt= '*g', label = 'v_y', markersize=5, capsize=5) - ax1.errorbar(modulation_depth_radial, 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='g') - h1, l1 = ax1.get_legend_handles_labels() - h2, l2 = ax2.get_legend_handles_labels() - ax1.legend(h1+h2, l1+l2, loc=0, prop={'size': 12, 'weight': 'bold'}) - 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(fx, fy, fz, dfx, dfy, dfz, v_x, v_y, v_z, modulation_depth, w_x, w_z, plot_against_mod_depth = True): - - w_xs = w_x * convert_modulation_depth_to_alpha(modulation_depth)[0] - new_aspect_ratio = w_xs / w_z - - 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(B_range = [0, 2.59]): - BField = np.arange(B_range[0], B_range[1], 1e-3) * u.G - a_s_array = np.zeros(len(BField)) * ac.a0 - for idx in range(len(BField)): - a_s_array[idx], a_bkg = scatteringLength(BField[idx]) - rmelmIdx = [i for i, x in enumerate(np.isinf(a_s_array.value)) if x] - for x in rmelmIdx: - a_s_array[x-1] = np.inf * ac.a0 - - plt.figure(figsize=(9, 7)) - plt.plot(BField, a_s_array/ac.a0, '-b') - plt.axhline(y = a_bkg/ac.a0, color = 'r', linestyle = '--') - plt.text(min(BField.value) + 0.5, (a_bkg/ac.a0).value + 1, '$a_{bkg}$ = %.2f a0' %((a_bkg/ac.a0).value), fontsize=14, fontweight='bold') - plt.xlim([min(BField.value), max(BField.value)]) - plt.ylim([65, 125]) - plt.xlabel('B field (G)', fontsize= 12, fontweight='bold') - plt.ylabel('Scattering length (a0)', fontsize= 12, fontweight='bold') - plt.tight_layout() - 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() - -##################################################################### \ No newline at end of file diff --git a/ODT Calculator/Helpers.py b/ODT Calculator/Helpers.py new file mode 100644 index 0000000..a7ee479 --- /dev/null +++ b/ODT Calculator/Helpers.py @@ -0,0 +1,80 @@ +import math +import numpy as np +from scipy import signal +from astropy import units as u, constants as ac + +DY_POLARIZABILITY = 184.4 # in a.u, most precise measured value of Dy polarizability +DY_MASS = 164*u.u +DY_DIPOLE_MOMENT = 9.93 * ac.muB + +##################################################################### +# HELPER FUNCTIONS # +##################################################################### + +def orderOfMagnitude(number): + return math.floor(math.log(number, 10)) + +def rotation_matrix(axis, theta): + """ + Return the rotation matrix associated with counterclockwise rotation about + the given axis by theta radians. + In 2-D it is just, + thetaInRadians = np.radians(theta) + c, s = np.cos(thetaInRadians), np.sin(thetaInRadians) + R = np.array(((c, -s), (s, c))) + In 3-D, one way to do it is use the Euler-Rodrigues Formula as is done here + """ + axis = np.asarray(axis) + axis = axis / math.sqrt(np.dot(axis, axis)) + a = math.cos(theta / 2.0) + b, c, d = -axis * math.sin(theta / 2.0) + aa, bb, cc, dd = a * a, b * b, c * c, d * d + bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d + return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)], + [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)], + [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) + +def find_nearest(array, value): + array = np.asarray(array) + idx = (np.abs(array - value)).argmin() + return idx + +def modulation_function(mod_amp, n_points, func = 'arccos'): + + if func == 'sin': + 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) + mod_func = mod_amp * np.concatenate((tmp_1, tmp_2)) + elif func == 'triangle': + phi = np.linspace(0, 2*np.pi, n_points) + mod_func = mod_amp * signal.sawtooth(phi, width = 0.5) # width of 0.5 gives symmetric rising triangle ramp + elif func == 'square': + phi = np.linspace(0, 1.99*np.pi, n_points) + mod_func = mod_amp * signal.square(phi, duty = 0.5) + else: + mod_func = None + + if mod_func is not None: + dx = (max(mod_func) - min(mod_func))/(2*n_points) + + return dx, mod_func + +##################################################################### +# BEAM PARAMETERS # +##################################################################### + +# Rayleigh length +def z_R(w_0, lamb)->np.ndarray: + return np.pi*w_0**2/lamb + +# Beam Radius +def w(pos, w_0, lamb): + return w_0*np.sqrt(1+(pos / z_R(w_0, lamb))**2) + + diff --git a/ODT Calculator/Plotter.py b/ODT Calculator/Plotter.py new file mode 100644 index 0000000..7e4a8c8 --- /dev/null +++ b/ODT Calculator/Plotter.py @@ -0,0 +1,342 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.ticker as mtick +from astropy import units as u, constants as ac + +##################################################################### +# 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 + happrox = harmonic_potential(Positions[axis, :].value, *popt) + fig = plt.figure(figsize=(12, 6)) + ax = fig.add_subplot(121) + ax.set_title('Fit to Potential') + plt.plot(Positions[axis, :].value, happrox, '-r', label = '\u03BD = %.1f \u00B1 %.2f Hz'% tuple([v,dv])) + 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([-TrapDepthsInKelvin[0].value, max(TrappingPotential[axis].value)]) + plt.grid(visible=1) + plt.legend(prop={'size': 12, 'weight': 'bold'}) + + bx = fig.add_subplot(122) + bx.set_title('Fit Residuals') + plt.plot(Positions[axis, :].value, TrappingPotential[axis].value - happrox, 'ob') + plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold') + plt.ylabel('$U_{trap} - U_{Harmonic}$', fontsize= 12, fontweight='bold') + plt.xlim([-10, 10]) + plt.ylim([-1e-2, 1e-2]) + plt.grid(visible=1) + plt.tight_layout() + plt.show() + +def plotGaussianFit(Positions, TrappingPotential, popt, pcov): + extracted_waist = popt[1] + dextracted_waist = pcov[1][1]**0.5 + gapprox = gaussian_potential(Positions, *popt) + fig = plt.figure(figsize=(12, 6)) + ax = fig.add_subplot(121) + ax.set_title('Fit to Potential') + plt.plot(Positions, gapprox, '-r', label = 'waist = %.1f \u00B1 %.2f um'% tuple([extracted_waist,dextracted_waist])) + plt.plot(Positions, TrappingPotential, 'ob', label = 'Gaussian Potential') + plt.xlabel('Distance (um)', fontsize= 12, fontweight='bold') + plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold') + plt.ylim([min(TrappingPotential), max(TrappingPotential)]) + plt.grid(visible=1) + plt.legend(prop={'size': 12, 'weight': 'bold'}) + + bx = fig.add_subplot(122) + bx.set_title('Fit Residuals') + plt.plot(Positions, TrappingPotential - gapprox, 'ob') + plt.xlabel('Distance (um)', 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 plotPotential(Positions, ComputedPotentials, options, Params = [], listToIterateOver = [], save = False): + + axis = options['axis'] + + plt.figure(figsize=(9, 7)) + for i in range(np.size(ComputedPotentials, 0)): + + if i % 2 == 0: + j = int(i / 2) + else: + 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] + + 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): + 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 - Horizontal' + elif axis == 1: + dir = 'Y - Propagation' + else: + dir = 'Z - Vertical' + + plt.ylim(top = 0) + plt.xlabel(dir + ' Direction (um)', fontsize= 12, fontweight='bold') + plt.ylabel('Trap Potential (uK)', fontsize= 12, fontweight='bold') + 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() + +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] + + 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() + +def plotAlphas(): + + 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) + + plt.figure() + plt.errorbar(fin_mod_dep, alpha_x, yerr = dalpha_x, fmt= 'ob', label = 'From Horz TF', markersize=5, capsize=5) + plt.errorbar(fin_mod_dep, alpha_y, yerr = dalpha_y, fmt= 'or', label = 'From Vert TF', 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.legend(prop={'size': 12, 'weight': 'bold'}) + plt.show() + +def plotTemperatures(w_x, w_z, plot_against_mod_depth = True): + + 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 + + plt.figure() + if plot_against_mod_depth: + plt.errorbar(fin_mod_dep, T_x, yerr = dT_x, fmt= 'ob', label = 'Horz direction', markersize=5, capsize=5) + plt.errorbar(fin_mod_dep, T_y, yerr = dT_y, fmt= 'or', label = 'Vert direction', 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', label = 'Horz direction', markersize=5, capsize=5) + plt.errorbar(measured_aspect_ratio, T_y, yerr = dT_y, fmt= 'or', label = 'Vert direction', markersize=5, capsize=5) + plt.plot(new_aspect_ratio, Temperatures, '--g') + xlabel = 'Aspect Ratio' + + plt.xlabel(xlabel, fontsize= 12, fontweight='bold') + plt.ylabel('Temperature (uK)', fontsize= 12, fontweight='bold') + plt.tight_layout() + plt.grid(visible=1) + plt.legend(prop={'size': 12, 'weight': 'bold'}) + 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, '-or', label = 'v_x') + ln2 = ax3.plot(modulation_depth, v_z, '-^b', label = 'v_z') + ax4 = ax3.twinx() + ln3 = ax4.plot(modulation_depth, v_y, '-*g', label = 'v_y') + xlabel = 'Modulation depth' + else: + ln1 = ax3.plot(new_aspect_ratio, v_x, '-or', 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, '-*g', 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='g') + 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(fx, dfx, fy, dfy, fz, dfz, modulation_depth_radial, modulation_depth_axial, w_x, w_z, plot_against_mod_depth = True): + + 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(modulation_depth_radial, fx, yerr = dfx, fmt= 'or', label = 'v_x', markersize=5, capsize=5) + ax2.errorbar(modulation_depth_axial, fy, yerr = dfy, fmt= '*g', label = 'v_y', markersize=5, capsize=5) + ax1.errorbar(modulation_depth_radial, 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='g') + h1, l1 = ax1.get_legend_handles_labels() + h2, l2 = ax2.get_legend_handles_labels() + ax1.legend(h1+h2, l1+l2, loc=0, prop={'size': 12, 'weight': 'bold'}) + 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(fx, fy, fz, dfx, dfy, dfz, v_x, v_y, v_z, modulation_depth, w_x, w_z, plot_against_mod_depth = True): + + w_xs = w_x * convert_modulation_depth_to_alpha(modulation_depth)[0] + new_aspect_ratio = w_xs / w_z + + 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(B_range = [0, 2.59]): + BField = np.arange(B_range[0], B_range[1], 1e-3) * u.G + a_s_array = np.zeros(len(BField)) * ac.a0 + for idx in range(len(BField)): + a_s_array[idx], a_bkg = scatteringLength(BField[idx]) + rmelmIdx = [i for i, x in enumerate(np.isinf(a_s_array.value)) if x] + for x in rmelmIdx: + a_s_array[x-1] = np.inf * ac.a0 + + plt.figure(figsize=(9, 7)) + plt.plot(BField, a_s_array/ac.a0, '-b') + plt.axhline(y = a_bkg/ac.a0, color = 'r', linestyle = '--') + plt.text(min(BField.value) + 0.5, (a_bkg/ac.a0).value + 1, '$a_{bkg}$ = %.2f a0' %((a_bkg/ac.a0).value), fontsize=14, fontweight='bold') + plt.xlim([min(BField.value), max(BField.value)]) + plt.ylim([65, 125]) + plt.xlabel('B field (G)', fontsize= 12, fontweight='bold') + plt.ylabel('Scattering length (a0)', fontsize= 12, fontweight='bold') + plt.tight_layout() + 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() + +##################################################################### \ No newline at end of file diff --git a/ODT Calculator/Potentials.py b/ODT Calculator/Potentials.py new file mode 100644 index 0000000..2bb5648 --- /dev/null +++ b/ODT Calculator/Potentials.py @@ -0,0 +1,105 @@ +import numpy as np +from astropy import units as u, constants as ac + +DY_POLARIZABILITY = 184.4 # in a.u, most precise measured value of Dy polarizability +DY_MASS = 164*u.u +DY_DIPOLE_MOMENT = 9.93 * ac.muB + +##################################################################### +# POTENTIALS # +##################################################################### + +def gravitational_potential(positions, m): + return m * ac.g0 * positions + +def single_gaussian_beam_potential(positions, waists, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um): + 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, waists, del_y, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um): + 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 modulated_single_gaussian_beam_potential(positions, waists, alpha = DY_POLARIZABILITY, P=1, wavelength=1.064*u.um, mod_amp=1): + mod_amp = mod_amp * waists[0] + n_points = len(positions[0,:]) + dx, x_mod = modulation_function(mod_amp, n_points, func = 'arccos') + 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) + dU = np.zeros(2*n_points) + for i in range(len(x_mod)): + dU = np.vstack((dU, np.exp(-2 * (np.subtract(x_mod[i], positions[0,:])/w(positions[1,:], waists[0], wavelength))**2))) + U = - U_tilde * A * 1/(2*mod_amp) * np.trapz(dU, dx = dx, axis = 0) + return U + +def harmonic_potential(pos, v, xoffset, yoffset, m = DY_MASS): + 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 + +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, waists, P, options, alpha = DY_POLARIZABILITY, wavelength=1.064*u.um): + + delta = options['delta'] + + foci_shift = options['foci_shift'] + focus_shift_beam_1 = foci_shift[0] + focus_shift_beam_2 = foci_shift[1] + + beam_disp = options['beam_disp'] + beam_1_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[0])).T * beam_disp[0].unit + beam_2_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[1])).T * beam_disp[1].unit + + beam_1_positions = positions + beam_1_disp + A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][0], wavelength)*w(beam_1_positions[1,:] + focus_shift_beam_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,:] + focus_shift_beam_1, waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + focus_shift_beam_1, waists[0][1], wavelength))**2)) + + R = rotation_matrix([0, 0, 1], np.radians(delta)) + beam_2_positions = np.dot(R, positions + beam_2_disp) + A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][0], wavelength)*w(beam_2_positions[1,:] + focus_shift_beam_2, 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,:] + focus_shift_beam_2, waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + focus_shift_beam_2, waists[1][1], wavelength))**2)) + + U = U_1 + U_2 + + return U + +def astigmatic_crossed_beam_potential(positions, waists, P, options, alpha = DY_POLARIZABILITY, wavelength=1.064*u.um): + + delta = options['delta'] + + del_y = options['foci_disp_crossed'] + del_y_1 = del_y[0] + del_y_2 = del_y[1] + + + foci_shift = options['foci_shift'] + focus_shift_beam_1 = foci_shift[0] + focus_shift_beam_2 = foci_shift[1] + + beam_disp = options['beam_disp'] + beam_1_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[0])).T * beam_disp[0].unit + beam_2_disp = (np.ones(np.shape(positions.T)) * np.array(beam_disp[1])).T * beam_disp[1].unit + + beam_1_positions = positions + beam_1_disp + A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] - (del_y_1/2) + focus_shift_beam_1, waists[0][0], wavelength)*w(beam_1_positions[1,:] + (del_y_1/2) + focus_shift_beam_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,:] - (del_y_1/2) + focus_shift_beam_1, waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + (del_y_1/2) + focus_shift_beam_1, waists[0][1], wavelength))**2)) + + R = rotation_matrix([0, 0, 1], np.radians(delta)) + beam_2_positions = np.dot(R, positions + beam_2_disp) + A_2 = 2*P[1]/(np.pi*w(beam_2_positions[1,:] - (del_y_2/2) + focus_shift_beam_2, waists[1][0], wavelength)*w(beam_2_positions[1,:] + (del_y_2/2) + focus_shift_beam_2, 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,:] - (del_y_2/2) + focus_shift_beam_2, waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + (del_y_2/2) + focus_shift_beam_2, waists[1][1], wavelength))**2)) + + U = U_1 + U_2 + + return U +