diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py index 42e2340..b88510f 100644 --- a/calculateDipoleTrapPotential.py +++ b/calculateDipoleTrapPotential.py @@ -1,6 +1,7 @@ import math import numpy as np import matplotlib.pyplot as plt +from scipy import signal from scipy.optimize import curve_fit from astropy import units as u, constants as ac @@ -37,15 +38,28 @@ def find_nearest(array, value): return idx def modulation_function(mod_amp, n_points, func = 'arccos'): - if func == 'arccos': + + if func == 'sin': phi = np.linspace(0, 2*np.pi, n_points) - first_half = 2/np.pi * np.arccos(phi/np.pi-1) - 1 - second_half = np.flip(first_half) - mod_func = mod_amp * np.concatenate((first_half, second_half)) - dx = (max(mod_func) - min(mod_func))/(2*n_points) - return dx, mod_func + mod_func = mod_amp * np.sin(phi) + elif func == 'arccos': + 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: - return None + 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 # @@ -146,13 +160,9 @@ def astigmatic_single_gaussian_beam_potential(positions: "np.ndarray|u.quantity. 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 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 - def modulated_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, mod_amp:"float|u.quantity.Quantity"=1)->np.ndarray: mod_amp = mod_amp * waists[0] - n_points = int(len(positions[0,:])/2) + 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) @@ -162,6 +172,14 @@ def modulated_single_gaussian_beam_potential(positions: "np.ndarray|u.quantity.Q 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 = 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 + +def gaussian_potential(pos, amp, waist, xoffset, yoffset): + U_Gaussian = amp * np.exp(-2 * ((pos + xoffset) / waist)**2) + yoffset + return U_Gaussian + ##################################################################### # COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS # ##################################################################### @@ -295,6 +313,23 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options): return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies +def extractWaist(Positions, TrappingPotential): + tmp_pos = Positions.value + tmp_pot = TrappingPotential.value + center_idx = np.argmin(tmp_pot) + + 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)) + xdata = tmp_pos[lb:ub] + Potential = tmp_pot[lb:ub] + + p0=[TrapMinimum, 30, TrapCenter, 0] + popt, pcov = curve_fit(gaussian_potential, xdata, Potential, p0) + return popt, pcov + ##################################################################### # PLOTTING # ##################################################################### @@ -303,15 +338,52 @@ def plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, axis, popt v = popt[0] dv = pcov[0][0]**0.5 happrox = harmonic_potential(Positions[axis, :].value, *popt) - plt.figure() + 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.tight_layout() 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_{Harmonic}$', 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): @@ -381,6 +453,7 @@ def plotIntensityProfileAndPotentials(Power, waists, alpha, wavelength, options) w_z = waists[1] extent = options['extent'] modulation = options['modulation'] + mod_func = options['modulation_function'] if not modulation: extent = 50 @@ -431,24 +504,36 @@ def plotIntensityProfileAndPotentials(Power, waists, alpha, wavelength, options) z_Positions = np.arange(-extent, extent, 1)*u.um mod_amp = mod_amp * w_x - n_points = int(len(x_Positions)/2) - dx, xmod_Positions = modulation_function(mod_amp, n_points, func = 'arccos') + 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') - A = 2*Power/(np.pi*w(0*u.um , w_x, wavelength)*w(0*u.um , w_z, wavelength)) + + ## 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 = np.transpose(intensity_profile[:, idx, :].to(u.MW/(u.cm*u.cm))) + 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 + 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) + fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(121) - ax.set_title('Intensity Profile ($MW/cm^2$)') - im = plt.imshow(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)]) + 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') @@ -456,8 +541,8 @@ def plotIntensityProfileAndPotentials(Power, waists, alpha, wavelength, options) bx = fig.add_subplot(122) bx.set_title('Trap Potential') - plt.plot(x_Positions, U[np.where(x_Positions==0)[0][0], :], label = 'X - Horizontal') - plt.plot(z_Positions, U[:, np.where(z_Positions==0)[0][0]], label = 'Z - Vertical') + 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') @@ -546,24 +631,30 @@ if __name__ == '__main__': """Plot transverse intensity profile and trap potential resulting for given parameters only""" # options = { - # 'extent': 70, # range of spatial coordinates in one direction to calculate trap potential over + # 'extent': 50, # range of spatial coordinates in one direction to calculate trap potential over # 'modulation': True, - # 'modulation_amplitude': 4.37 + # 'modulation_function': 'arccos', + # 'modulation_amplitude': 2.12 # } # plotIntensityProfileAndPotentials(Power, [w_x, w_z], Polarizability, Wavelength, options) + """Plot gaussian fit for trap potential resulting from modulation for given parameters only""" + # plotGaussianFit(x_Positions, x_Potential, poptx, pcovx) + # plotGaussianFit(z_Positions, z_Potential, poptx, pcovx) """Calculate relevant parameters for evaporative cooling""" - AtomNumber = 1.13 * 1e7 - Temperature = 100 * u.uK - BField = 2.1 * u.G - modulation = False + AtomNumber = 1.00 * 1e7 + BField = 2.5 * u.G + modulation = True 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) @@ -585,4 +676,5 @@ if __name__ == '__main__': # plotScatteringLengths() + \ No newline at end of file