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