karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Added functionality to extract aspect ratio through a Gaussian fit, added plotting of fit residuals and other small bugfixes, changes

Browse Source
This commit is contained in:
Karthik 2023-02-17 19:23:53 +01:00
parent e7d2255ac9
commit 64132018e2
1 changed files with 120 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
calculateDipoleTrapPotential.py
View File

@ -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()