Added astigmatic crossed dipole trap, corrected inconsistency in the naming and usage of the variable for polarizability.

This commit is contained in:
Karthik 2023-03-14 19:26:43 +01:00
parent 2fdf592dda
commit 6f018c0281
2 changed files with 270 additions and 411 deletions

File diff suppressed because one or more lines are too long

View File

@ -87,10 +87,10 @@ def calculateAtomNumber(NCount, pixel_size = 5.86 * u.um, magnification = 0.5, e
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')
def particleDensity(w_x, w_z, Power, N, T, m = 164*u.u): # For a thermal cloud
v_x = calculateTrapFrequency(w_x, w_z, Power, dir = 'x')
v_y = calculateTrapFrequency(w_x, w_z, Power, dir = 'y')
v_z = calculateTrapFrequency(w_x, w_z, Power, 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 calculateParticleDensityFromMeasurements(v_x, dv_x, v_y, dv_y, v_z, dv_z, w_x, w_z, T_x, T_y, dT_x, dT_y, modulation_depth, N, m = 164*u.u):
@ -167,11 +167,11 @@ def dipolarLength(mu = 9.93 * ac.muB, m = 164*u.u):
def scatteringCrossSection(B):
return 8 * np.pi * scatteringLength(B)[0]**2 + ((32*np.pi)/45) * dipolarLength()**2
def calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N, T, B): #For a 3D Harmonic Trap
return (particleDensity(w_x, w_z, Power, Polarizability, N, T) * scatteringCrossSection(B) * meanThermalVelocity(T) / (2 * np.sqrt(2))).decompose()
def calculateElasticCollisionRate(w_x, w_z, Power, N, T, B): #For a 3D Harmonic Trap
return (particleDensity(w_x, w_z, Power, N, T) * scatteringCrossSection(B) * meanThermalVelocity(T) / (2 * np.sqrt(2))).decompose()
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 calculatePSD(w_x, w_z, Power, N, T):
return (particleDensity(w_x, w_z, Power, N, T) * 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]
@ -213,19 +213,19 @@ def convert_modulation_depth_to_temperature(modulation_depth):
def gravitational_potential(positions, m):
return m * ac.g0 * positions
def single_gaussian_beam_potential(positions, waists, alpha, P=1, wavelength=1.064*u.um):
def single_gaussian_beam_potential(positions, waists, alpha = 184.4, 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, P=1, wavelength=1.064*u.um):
def astigmatic_single_gaussian_beam_potential(positions, waists, del_y, alpha = 184.4, 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, P=1, wavelength=1.064*u.um, mod_amp=1):
def modulated_single_gaussian_beam_potential(positions, waists, alpha = 184.4, 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')
@ -245,7 +245,7 @@ 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):
def crossed_beam_potential(positions, theta, waists, P, alpha = 184.4, 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))
@ -262,15 +262,35 @@ def crossed_beam_potential(positions, theta, waists, P, alpha, wavelength=1.064*
return U
def astigmatic_crossed_beam_potential(positions, theta, waists, P, del_y, alpha = 184.4, wavelength=1.064*u.um):
del_y_1 = del_y[0]
del_y_2 = del_y[1]
beam_1_positions = positions
A_1 = 2*P[0]/(np.pi*w(beam_1_positions[1,:] - (del_y_1/2), waists[0][0], wavelength)*w(beam_1_positions[1,:] + (del_y_1/2), 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), waists[0][0], wavelength))**2 + (beam_1_positions[2,:]/w(beam_1_positions[1,:] + (del_y_1/2), 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,:] - (del_y_2/2), waists[1][0], wavelength)*w(beam_2_positions[1,:] + (del_y_2/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), waists[1][0], wavelength))**2 + (beam_2_positions[2,:]/w(beam_2_positions[1,:] + (del_y_2/2), waists[1][1], wavelength))**2))
U = U_1 + U_2
return U
#####################################################################
# COMPUTE/EXTRACT TRAP POTENTIAL AND PARAMETERS #
#####################################################################
def trap_depth(w_1, w_2, P, alpha):
def trap_depth(w_1, w_2, P, alpha = 184.4):
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, dir = 'x', m = 164*u.u):
TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability)
def calculateTrapFrequency(w_x, w_z, Power, dir = 'x', m = 164*u.u):
TrapDepth = trap_depth(w_x, w_z, Power)
TrapFrequency = np.nan
if dir == 'x':
TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth / (m*w_x**2))).decompose()
@ -281,10 +301,10 @@ def calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x', m = 164*u
TrapFrequency = ((1/(2 * np.pi)) * np.sqrt(4 * TrapDepth/ (m*w_z**2))).decompose()
return round(TrapFrequency.value, 2)*u.Hz
def calculateCrossedBeamTrapFrequency(delta, Waists, Powers, dir = 'x', m = 164*u.u, Polarizability = 184.4, wavelength=1.064*u.um):
def calculateCrossedBeamTrapFrequency(delta, Waists, Powers, dir = 'x', m = 164*u.u, wavelength=1.064*u.um):
TrapDepth_1 = trap_depth(Waists[0][0], Waists[1][0], Powers[0], alpha=Polarizability)
TrapDepth_2 = trap_depth(Waists[0][1], Waists[1][1], Powers[1], alpha=Polarizability)
TrapDepth_1 = trap_depth(Waists[0][0], Waists[1][0], Powers[0])
TrapDepth_2 = trap_depth(Waists[0][1], Waists[1][1], Powers[1])
w_x1 = Waists[0][0]
w_z1 = Waists[1][0]
@ -292,8 +312,8 @@ def calculateCrossedBeamTrapFrequency(delta, Waists, Powers, dir = 'x', m = 164*
w_y2 = Waists[0][1]
w_z2 = Waists[1][1]
zReff_1 = np.sqrt(2) * z_R(w_x1, 1.064*u.um) * z_R(w_z1, 1.064*u.um) / np.sqrt(z_R(w_x1, 1.064*u.um)**2 + z_R(w_z1, 1.064*u.um)**2)
zReff_2 = np.sqrt(2) * z_R(w_y2, 1.064*u.um) * z_R(w_z2, 1.064*u.um) / np.sqrt(z_R(w_y2, 1.064*u.um)**2 + z_R(w_z2, 1.064*u.um)**2)
zReff_1 = np.sqrt(2) * z_R(w_x1, wavelength) * z_R(w_z1, wavelength) / np.sqrt(z_R(w_x1, wavelength)**2 + z_R(w_z1, wavelength)**2)
zReff_2 = np.sqrt(2) * z_R(w_y2, wavelength) * z_R(w_z2, wavelength) / np.sqrt(z_R(w_y2, wavelength)**2 + z_R(w_z2, wavelength)**2)
wy2alpha = np.sqrt(((np.cos(np.radians(90 - delta)) / w_y2)**2 + (np.sin(np.radians(90 - delta)) / (2 * zReff_2))**2)**(-1))
zR2alpha = np.sqrt(((np.sin(np.radians(90 - delta)) / w_y2)**2 + (np.cos(np.radians(90 - delta)) / (2 * zReff_2))**2)**(-1))
@ -328,7 +348,7 @@ def extractTrapFrequency(Positions, TrappingPotential, axis):
dv = pcov[0][0]**0.5
return v, dv, popt, pcov
def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
def computeTrapPotential(w_x, w_z, Power, options):
axis = options['axis']
extent = options['extent']
@ -343,7 +363,7 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
w_x = w_x * (aspect_ratio / current_ar)
TrappingPotential = []
TrapDepth = trap_depth(w_x, w_z, Power, alpha=Polarizability)
TrapDepth = trap_depth(w_x, w_z, Power)
IdealTrapDepthInKelvin = (TrapDepth/ac.k_B).to(u.uK)
projection_axis = np.array([0, 1, 0]) # default
@ -366,74 +386,68 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
Positions = np.vstack((x_Positions, y_Positions, z_Positions)) * projection_axis[:, np.newaxis]
if not crossed:
IdealTrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability)
IdealTrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power)
IdealTrappingPotential = IdealTrappingPotential * (np.ones((3, len(IdealTrappingPotential))) * projection_axis[:, np.newaxis])
IdealTrappingPotential = (IdealTrappingPotential/ac.k_B).to(u.uK)
else:
delta = options['delta']
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, delta, 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
m = 164*u.u
gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity']
theta = options['theta']
tilt_axis = options['tilt_axis']
if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, alpha = Polarizability)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
if gravity and not astigmatism:
# Influence of Gravity
m = 164*u.u
gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity']
theta = options['theta']
tilt_axis = options['tilt_axis']
if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif not gravity and astigmatism:
# Influence of Astigmatism
disp_foci = options['disp_foci']
TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis])
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif not gravity and astigmatism:
# Influence of Astigmatism
disp_foci = options['disp_foci']
TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci, alpha = Polarizability)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis])
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif gravity and astigmatism:
# Influence of Gravity and Astigmatism
m = 164*u.u
gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity']
theta = options['theta']
tilt_axis = options['tilt_axis']
disp_foci = options['disp_foci']
if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci, alpha = Polarizability)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
else:
TrappingPotential = IdealTrappingPotential
if not crossed:
elif gravity and astigmatism:
# Influence of Gravity and Astigmatism
m = 164*u.u
gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity']
theta = options['theta']
tilt_axis = options['tilt_axis']
disp_foci = options['disp_foci']
if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = astigmatic_single_gaussian_beam_potential(Positions, np.asarray([w_x.value, w_z.value])*u.um, P = Power, del_y = disp_foci)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
else:
TrappingPotential = IdealTrappingPotential
infls = np.where(np.diff(np.sign(np.gradient(np.gradient(TrappingPotential[axis].value)))))[0]
if TrappingPotential[axis][0] > TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][infls[1]:-1]) - min(TrappingPotential[axis][infls[0]:infls[1]])
elif TrappingPotential[axis][0] < TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][0:infls[0]]) - min(TrappingPotential[axis][infls[0]:infls[1]])
else:
EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin
try:
if TrappingPotential[axis][0] > TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][infls[1]:-1]) - min(TrappingPotential[axis][infls[0]:infls[1]])
elif TrappingPotential[axis][0] < TrappingPotential[axis][-1]:
EffectiveTrapDepthInKelvin = max(TrappingPotential[axis][0:infls[0]]) - min(TrappingPotential[axis][infls[0]:infls[1]])
else:
EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin
except:
EffectiveTrapDepthInKelvin = np.nan
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')
v_x = calculateTrapFrequency(w_x, w_z, Power, dir = 'x')
v_y = calculateTrapFrequency(w_x, w_z, Power, dir = 'y')
v_z = calculateTrapFrequency(w_x, w_z, Power, dir = 'z')
CalculatedTrapFrequencies = [v_x, v_y, v_z]
v, dv, popt, pcov = extractTrapFrequency(Positions, IdealTrappingPotential, axis)
@ -458,6 +472,53 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies
else:
delta = options['delta']
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, delta, waists, P = Power)
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
m = 164*u.u
gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity']
theta = options['theta']
tilt_axis = options['tilt_axis']
if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = crossed_beam_potential(Positions, delta, waists, P = Power)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif not gravity and astigmatism:
# Influence of Astigmatism
disp_foci = options['disp_foci_crossed']
TrappingPotential = astigmatic_crossed_beam_potential(Positions, delta, waists, P = Power, del_y = disp_foci)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis])
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
elif gravity and astigmatism:
# Influence of Gravity and Astigmatism
m = 164*u.u
gravity_axis = np.array([0, 0, 1])
tilt_gravity = options['tilt_gravity']
theta = options['theta']
tilt_axis = options['tilt_axis']
disp_foci = options['disp_foci_crossed']
if tilt_gravity:
R = rotation_matrix(tilt_axis, np.radians(theta))
gravity_axis = np.dot(R, gravity_axis)
gravity_axis_positions = np.vstack((x_Positions, y_Positions, z_Positions)) * gravity_axis[:, np.newaxis]
TrappingPotential = astigmatic_crossed_beam_potential(Positions, delta, waists, P = Power, del_y = disp_foci)
TrappingPotential = TrappingPotential * (np.ones((3, len(TrappingPotential))) * projection_axis[:, np.newaxis]) + gravitational_potential(gravity_axis_positions, m)
TrappingPotential = (TrappingPotential/ac.k_B).to(u.uK)
else:
TrappingPotential = IdealTrappingPotential
return Positions, TrappingPotential
def extractWaist(Positions, TrappingPotential):
@ -477,7 +538,7 @@ def extractWaist(Positions, TrappingPotential):
popt, pcov = curve_fit(gaussian_potential, xdata, Potential, p0)
return popt, pcov
def computeIntensityProfileAndPotentials(Power, waists, alpha, wavelength, options):
def computeIntensityProfileAndPotentials(Power, waists, wavelength, options, alpha = 184.4):
w_x = waists[0]
w_z = waists[1]
extent = options['extent']
@ -492,7 +553,6 @@ def computeIntensityProfileAndPotentials(Power, waists, alpha, wavelength, optio
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')
@ -877,4 +937,6 @@ def plotCollisionRatesAndPSD(Gamma_elastic, PSD, modulation_depth, new_aspect_ra
plt.grid(visible=1)
plt.show()
#####################################################################
#####################################################################
#Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability