Script to estimate trap frequencies of the hybrid trap of AL, cODT.
This commit is contained in:
parent
a3899a7827
commit
750bd1bf69
368
Estimations/EstimatesForHybridTrap.m
Normal file
368
Estimations/EstimatesForHybridTrap.m
Normal file
@ -0,0 +1,368 @@
|
||||
%% Physical constants
|
||||
PlanckConstant = 6.62607015E-34;
|
||||
PlanckConstantReduced = 6.62607015E-34/(2*pi);
|
||||
FineStructureConstant = 7.2973525698E-3;
|
||||
ElectronMass = 9.10938291E-31;
|
||||
GravitationalConstant = 6.67384E-11;
|
||||
ProtonMass = 1.672621777E-27;
|
||||
AtomicMassUnit = 1.660539066E-27;
|
||||
BohrRadius = 5.2917721067E-11;
|
||||
BohrMagneton = 9.274009994E-24;
|
||||
BoltzmannConstant = 1.38064852E-23;
|
||||
StandardGravityAcceleration = 9.80665;
|
||||
SpeedOfLight = 299792458;
|
||||
StefanBoltzmannConstant = 5.670373E-8;
|
||||
ElectronCharge = 1.602176634E-19;
|
||||
VacuumPermeability = 1.25663706212E-6;
|
||||
DielectricConstant = 8.8541878128E-12;
|
||||
ElectronGyromagneticFactor = -2.00231930436153;
|
||||
AvogadroConstant = 6.02214076E23;
|
||||
ZeroKelvin = 273.15;
|
||||
GravitationalAcceleration = 9.80553;
|
||||
VacuumPermittivity = 1 / (SpeedOfLight^2 * VacuumPermeability);
|
||||
HartreeEnergy = ElectronCharge^2 / (4 * pi * VacuumPermittivity * BohrRadius);
|
||||
AtomicUnitOfPolarizability = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3
|
||||
|
||||
% Dy specific constants
|
||||
Dy164Mass = 163.929174751*AtomicMassUnit;
|
||||
Dy164IsotopicAbundance = 0.2826;
|
||||
DyMagneticMoment = 9.93*BohrMagneton;
|
||||
|
||||
%% Helper Functions
|
||||
|
||||
function R = rotation_matrix(axis, theta)
|
||||
% rotation_matrix returns the 3x3 rotation matrix associated with
|
||||
% counterclockwise rotation about the given axis by theta radians.
|
||||
|
||||
% Normalize the axis
|
||||
axis = axis / sqrt(dot(axis, axis));
|
||||
|
||||
% Compute the rotation matrix components
|
||||
a = cos(theta / 2.0);
|
||||
b = -axis(1) * sin(theta / 2.0);
|
||||
c = -axis(2) * sin(theta / 2.0);
|
||||
d = -axis(3) * sin(theta / 2.0);
|
||||
|
||||
% Precompute terms for readability
|
||||
aa = a * a; bb = b * b; cc = c * c; dd = d * d;
|
||||
bc = b * c; ad = a * d; ac = a * c; ab = a * b;
|
||||
bd = b * d; cd = c * d;
|
||||
|
||||
% Construct the rotation matrix
|
||||
R = [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];
|
||||
end
|
||||
|
||||
%% Potentials
|
||||
|
||||
function U_gravity = gravitational_potential(positions, m)
|
||||
g = 9.81; % Earth's gravitational constant
|
||||
U = m * g * positions;
|
||||
end
|
||||
|
||||
function U_Harmonic = harmonic_potential(pos, v, xoffset, yoffset, m)
|
||||
% Calculate the harmonic potential
|
||||
k_B = physical_constant('Boltzmann constant');
|
||||
Hz = 1; % for unit Hz
|
||||
U_Harmonic = ((0.5 * m * (2 * pi * v * Hz).^2 * (pos - xoffset).^2) / k_B) + yoffset;
|
||||
end
|
||||
|
||||
function U_crossedodt = crossed_beam_potential(positions, waists, P, options, alpha)
|
||||
% Crossed beam potential
|
||||
if nargin < 5
|
||||
alpha = 184.4; % Dy polarizability default
|
||||
end
|
||||
wavelength = 1.064e-6; % wavelength in meters
|
||||
|
||||
delta = options.delta;
|
||||
foci_shift = options.foci_shift;
|
||||
focus_shift_beam_1 = foci_shift(1);
|
||||
focus_shift_beam_2 = foci_shift(2);
|
||||
|
||||
beam_disp = options.beam_disp;
|
||||
beam_1_disp = beam_disp(1);
|
||||
beam_2_disp = beam_disp(2);
|
||||
|
||||
beam_1_positions = positions + beam_1_disp;
|
||||
A_1 = 2 * P(1) ./ (pi * w(beam_1_positions(2,:) + focus_shift_beam_1, waists(1,1), wavelength) ...
|
||||
.* w(beam_1_positions(2,:) + focus_shift_beam_1, waists(1,2), wavelength));
|
||||
U_1_tilde = (1 / (2 * 8.854e-12 * 3e8)) * alpha * (4 * pi * 8.854e-12 * 5.29177e-11^3);
|
||||
U_1 = -U_1_tilde * A_1 .* exp(-2 * ((beam_1_positions(1,:) ./ w(beam_1_positions(2,:) + focus_shift_beam_1, waists(1,1), wavelength)).^2 ...
|
||||
+ (beam_1_positions(3,:) ./ w(beam_1_positions(2,:) + focus_shift_beam_1, waists(1,2), wavelength)).^2));
|
||||
|
||||
R = rotation_matrix([0, 0, 1], deg2rad(delta));
|
||||
beam_2_positions = R * (positions + beam_2_disp);
|
||||
A_2 = 2 * P(2) ./ (pi * w(beam_2_positions(2,:) + focus_shift_beam_2, waists(2,1), wavelength) ...
|
||||
.* w(beam_2_positions(2,:) + focus_shift_beam_2, waists(2,2), wavelength));
|
||||
U_2_tilde = (1 / (2 * 8.854e-12 * 3e8)) * alpha * (4 * pi * 8.854e-12 * 5.29177e-11^3);
|
||||
U_2 = -U_2_tilde * A_2 .* exp(-2 * ((beam_2_positions(1,:) ./ w(beam_2_positions(2,:) + focus_shift_beam_2, waists(2,1), wavelength)).^2 ...
|
||||
+ (beam_2_positions(3,:) ./ w(beam_2_positions(2,:) + focus_shift_beam_2, waists(2,2), wavelength)).^2));
|
||||
|
||||
U = U_1 + U_2;
|
||||
end
|
||||
|
||||
function U_astigcrossedodt = astigmatic_crossed_beam_potential(positions, waists, P, options, alpha, wavelength)
|
||||
% ASTIGMATIC_CROSSED_BEAM_POTENTIAL - Calculate the potential of an astigmatic crossed beam trap.
|
||||
%
|
||||
% U = ASTIGMATIC_CROSSED_BEAM_POTENTIAL(positions, waists, P, options, alpha, wavelength)
|
||||
%
|
||||
% Parameters:
|
||||
% positions - 3xN matrix of positions for the particle/atom.
|
||||
% waists - 2x2 cell array of beam waists, where each element is a 2-element vector [waist_x, waist_y].
|
||||
% P - 1x2 array of powers for the two beams.
|
||||
% options - Struct containing options:
|
||||
% - delta: Angle between beams in degrees.
|
||||
% - foci_disp_crossed: Displacement in the y-direction for both beams [del_y_1, del_y_2].
|
||||
% - foci_shift: Foci shifts for both beams [shift_1, shift_2].
|
||||
% - beam_disp: Displacement for both beams.
|
||||
% alpha - Polarizability of the atom/particle.
|
||||
% wavelength- Wavelength of the beams.
|
||||
%
|
||||
% Returns:
|
||||
% U - The total potential at the given positions.
|
||||
|
||||
if nargin < 5
|
||||
alpha = DY_POLARIZABILITY;
|
||||
end
|
||||
if nargin < 6
|
||||
wavelength = 1.064 * 1e-6; % Default wavelength in meters (1.064 um)
|
||||
end
|
||||
|
||||
delta = options.delta;
|
||||
|
||||
del_y = options.foci_disp_crossed;
|
||||
del_y_1 = del_y(1);
|
||||
del_y_2 = del_y(2);
|
||||
|
||||
foci_shift = options.foci_shift;
|
||||
focus_shift_beam_1 = foci_shift(1);
|
||||
focus_shift_beam_2 = foci_shift(2);
|
||||
|
||||
beam_disp = options.beam_disp;
|
||||
beam_1_disp = repmat(beam_disp{1}, 1, size(positions, 2));
|
||||
beam_2_disp = repmat(beam_disp{2}, 1, size(positions, 2));
|
||||
|
||||
% Calculate beam 1 potential
|
||||
beam_1_positions = positions + beam_1_disp;
|
||||
A_1 = 2*P(1) / (pi * w(beam_1_positions(2,:) - (del_y_1/2) + focus_shift_beam_1, waists{1}(1), wavelength) ...
|
||||
* w(beam_1_positions(2,:) + (del_y_1/2) + focus_shift_beam_1, waists{1}(2), wavelength));
|
||||
U_1_tilde = (1 / (2 * eps0 * c)) * alpha * (4 * pi * eps0 * a0^3);
|
||||
U_1 = - U_1_tilde * A_1 .* exp(-2 * ((beam_1_positions(1,:) ./ w(beam_1_positions(2,:) - (del_y_1/2) + focus_shift_beam_1, waists{1}(1), wavelength)).^2 ...
|
||||
+ (beam_1_positions(3,:) ./ w(beam_1_positions(2,:) + (del_y_1/2) + focus_shift_beam_1, waists{1}(2), wavelength)).^2));
|
||||
|
||||
% Rotation matrix for beam 2
|
||||
R = rotation_matrix([0, 0, 1], deg2rad(delta));
|
||||
beam_2_positions = R * (positions + beam_2_disp);
|
||||
|
||||
% Calculate beam 2 potential
|
||||
A_2 = 2*P(2) / (pi * w(beam_2_positions(2,:) - (del_y_2/2) + focus_shift_beam_2, waists{2}(1), wavelength) ...
|
||||
* w(beam_2_positions(2,:) + (del_y_2/2) + focus_shift_beam_2, waists{2}(2), wavelength));
|
||||
U_2_tilde = (1 / (2 * eps0 * c)) * alpha * (4 * pi * eps0 * a0^3);
|
||||
U_2 = - U_2_tilde * A_2 .* exp(-2 * ((beam_2_positions(1,:) ./ w(beam_2_positions(2,:) - (del_y_2/2) + focus_shift_beam_2, waists{2}(1), wavelength)).^2 ...
|
||||
+ (beam_2_positions(3,:) ./ w(beam_2_positions(2,:) + (del_y_2/2) + focus_shift_beam_2, waists{2}(2), wavelength)).^2));
|
||||
|
||||
% Total potential
|
||||
U = U_1 + U_2;
|
||||
end
|
||||
|
||||
%% Estimators
|
||||
|
||||
function depth = trap_depth(w_1, w_2, P, alpha)
|
||||
% Constants
|
||||
if nargin < 4
|
||||
alpha = DY_POLARIZABILITY; % Default polarizability if not provided
|
||||
end
|
||||
|
||||
% Define constants
|
||||
eps0 = 8.854187817e-12; % Permittivity of free space (F/m)
|
||||
c = 299792458; % Speed of light in vacuum (m/s)
|
||||
a0 = 5.29177210903e-11; % Bohr radius (m)
|
||||
|
||||
% Calculate the trap depth
|
||||
depth = (2 * P) / (pi * w_1 * w_2) * (1 / (2 * eps0 * c)) * alpha * (4 * pi * eps0 * a0^3);
|
||||
end
|
||||
|
||||
function [v, dv, popt, pcov] = extractTrapFrequency(Positions, TrappingPotential, axis)
|
||||
% EXTRACTTRAPFREQUENCY - Extract the trapping frequency by fitting a harmonic potential.
|
||||
%
|
||||
% [v, dv, popt, pcov] = EXTRACTTRAPFREQUENCY(Positions, TrappingPotential, axis)
|
||||
%
|
||||
% Parameters:
|
||||
% Positions - 2D matrix containing the positions of the atoms/particles.
|
||||
% TrappingPotential - Vector containing the potential values for the corresponding positions.
|
||||
% axis - Index of the axis (1, 2, or 3) to extract the frequency.
|
||||
%
|
||||
% Returns:
|
||||
% v - Extracted trapping frequency.
|
||||
% dv - Error/uncertainty in the frequency.
|
||||
% popt - Optimized parameters from the curve fitting.
|
||||
% pcov - Covariance matrix from the curve fitting.
|
||||
|
||||
tmp_pos = Positions(axis, :);
|
||||
tmp_pot = TrappingPotential(axis, :);
|
||||
|
||||
% Find the index of the potential's minimum (trap center)
|
||||
[~, center_idx] = min(tmp_pot);
|
||||
|
||||
% Define a range around the center for fitting (1/100th of the total length)
|
||||
lb = round(center_idx - length(tmp_pot)/200);
|
||||
ub = round(center_idx + length(tmp_pot)/200);
|
||||
|
||||
% Ensure boundaries are valid
|
||||
lb = max(lb, 1); % Lower bound cannot be less than 1
|
||||
ub = min(ub, length(tmp_pot)); % Upper bound cannot exceed length
|
||||
|
||||
% Extract the subset of position and potential data around the center
|
||||
xdata = tmp_pos(lb:ub);
|
||||
Potential = tmp_pot(lb:ub);
|
||||
|
||||
% Initial guess for fitting parameters [v, x0, amplitude]
|
||||
p0 = [1e3, tmp_pos(center_idx), -100];
|
||||
|
||||
% Perform curve fitting using harmonic potential
|
||||
try
|
||||
% Define harmonic potential function (anonymous function)
|
||||
harmonic_potential = @(p, x) p(3) + 0.5 * p(1) * (x - p(2)).^2;
|
||||
|
||||
% Fit the potential data to the harmonic potential model
|
||||
[popt, ~, pcov] = lsqcurvefit(harmonic_potential, p0, xdata, Potential);
|
||||
|
||||
v = popt(1); % Extract fitted frequency
|
||||
dv = sqrt(pcov(1, 1)); % Uncertainty in the frequency
|
||||
catch
|
||||
% In case of fitting failure, return NaNs
|
||||
popt = NaN;
|
||||
pcov = NaN;
|
||||
v = NaN;
|
||||
dv = NaN;
|
||||
end
|
||||
end
|
||||
|
||||
function [Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, ExtractedTrapFrequencies] = computeTrapPotential(w_x, w_z, Power, options)
|
||||
|
||||
% Unpack the options structure
|
||||
axis = options.axis;
|
||||
extent = options.extent;
|
||||
gravity = options.gravity;
|
||||
astigmatism = options.astigmatism;
|
||||
modulation = options.modulation;
|
||||
|
||||
% Apply modulation if necessary
|
||||
if modulation
|
||||
aspect_ratio = options.aspect_ratio;
|
||||
current_ar = w_x / w_z;
|
||||
w_x = w_x * (aspect_ratio / current_ar);
|
||||
end
|
||||
|
||||
% Compute ideal trap depth
|
||||
TrapDepth = trap_depth(w_x, w_z, Power);
|
||||
IdealTrapDepthInKelvin = TrapDepth / ac.k_B * u.uK; % Conversion to uK
|
||||
|
||||
% Default projection axis
|
||||
projection_axis = [0; 1; 0]; % Default Y-axis
|
||||
|
||||
% Set projection axis based on the given axis option
|
||||
if axis == 0
|
||||
projection_axis = [1; 0; 0]; % X-axis
|
||||
elseif axis == 1
|
||||
projection_axis = [0; 1; 0]; % Y-axis
|
||||
elseif axis == 2
|
||||
projection_axis = [0; 0; 1]; % Z-axis
|
||||
else
|
||||
projection_axis = [1; 1; 1]; % Diagonal
|
||||
end
|
||||
|
||||
% Define positions along the specified extent
|
||||
x_Positions = (-extent:0.05:extent) * u.um;
|
||||
y_Positions = (-extent:0.05:extent) * u.um;
|
||||
z_Positions = (-extent:0.05:extent) * u.um;
|
||||
|
||||
% Stack positions along the projection axis
|
||||
Positions = [x_Positions; y_Positions; z_Positions] .* projection_axis;
|
||||
|
||||
% Define the waists for the beams
|
||||
waists = [w_x(1), w_z(1); w_x(2), w_z(2)] * u.um;
|
||||
|
||||
% Calculate the ideal trapping potential
|
||||
IdealTrappingPotential = crossed_beam_potential(Positions, waists, Power, options);
|
||||
IdealTrappingPotential = IdealTrappingPotential .* (ones(3, length(IdealTrappingPotential)) .* projection_axis);
|
||||
IdealTrappingPotential = IdealTrappingPotential / ac.k_B * u.uK; % Convert to uK
|
||||
|
||||
% Initialize the TrappingPotential variable
|
||||
TrappingPotential = [];
|
||||
|
||||
% Handle different combinations of gravity and astigmatism
|
||||
if gravity && ~astigmatism
|
||||
% Gravity effect only
|
||||
m = DY_MASS;
|
||||
gravity_axis = [0; 0; 1]; % Z-axis for gravity
|
||||
if options.tilt_gravity
|
||||
R = rotation_matrix(options.tilt_axis, deg2rad(options.theta));
|
||||
gravity_axis = R * gravity_axis;
|
||||
end
|
||||
gravity_axis_positions = [x_Positions; y_Positions; z_Positions] .* gravity_axis;
|
||||
TrappingPotential = crossed_beam_potential(Positions, waists, Power, options);
|
||||
TrappingPotential = TrappingPotential .* (ones(3, length(TrappingPotential)) .* projection_axis) + gravitational_potential(gravity_axis_positions, m);
|
||||
TrappingPotential = TrappingPotential / ac.k_B * u.uK; % Convert to uK
|
||||
|
||||
elseif ~gravity && astigmatism
|
||||
% Astigmatism effect only
|
||||
TrappingPotential = astigmatic_crossed_beam_potential(Positions, waists, Power, options);
|
||||
TrappingPotential = TrappingPotential .* (ones(3, length(TrappingPotential)) .* projection_axis);
|
||||
TrappingPotential = TrappingPotential / ac.k_B * u.uK; % Convert to uK
|
||||
|
||||
elseif gravity && astigmatism
|
||||
% Both gravity and astigmatism
|
||||
m = DY_MASS;
|
||||
gravity_axis = [0; 0; 1]; % Z-axis for gravity
|
||||
if options.tilt_gravity
|
||||
R = rotation_matrix(options.tilt_axis, deg2rad(options.theta));
|
||||
gravity_axis = R * gravity_axis;
|
||||
end
|
||||
gravity_axis_positions = [x_Positions; y_Positions; z_Positions] .* gravity_axis;
|
||||
TrappingPotential = astigmatic_crossed_beam_potential(Positions, waists, Power, options);
|
||||
TrappingPotential = TrappingPotential .* (ones(3, length(TrappingPotential)) .* projection_axis) + gravitational_potential(gravity_axis_positions, m);
|
||||
TrappingPotential = TrappingPotential / ac.k_B * u.uK; % Convert to uK
|
||||
|
||||
else
|
||||
% Default to ideal trapping potential
|
||||
TrappingPotential = IdealTrappingPotential;
|
||||
end
|
||||
|
||||
% Calculate inflection points and effective trap depth
|
||||
infls = find(diff(sign(gradient(gradient(TrappingPotential(axis,:))))));
|
||||
|
||||
try
|
||||
if TrappingPotential(axis, 1) > TrappingPotential(axis, end)
|
||||
EffectiveTrapDepthInKelvin = max(TrappingPotential(axis, infls(2):end)) - min(TrappingPotential(axis, infls(1):infls(2)));
|
||||
elseif TrappingPotential(axis, 1) < TrappingPotential(axis, end)
|
||||
EffectiveTrapDepthInKelvin = max(TrappingPotential(axis, 1:infls(1))) - min(TrappingPotential(axis, infls(1):infls(2)));
|
||||
else
|
||||
EffectiveTrapDepthInKelvin = IdealTrapDepthInKelvin;
|
||||
end
|
||||
catch
|
||||
EffectiveTrapDepthInKelvin = NaN;
|
||||
end
|
||||
|
||||
% Trap depths in Kelvin
|
||||
TrapDepthsInKelvin = [IdealTrapDepthInKelvin, EffectiveTrapDepthInKelvin];
|
||||
|
||||
% Extract trap frequencies from the ideal trapping potential
|
||||
[v, dv, popt, pcov] = extractTrapFrequency(Positions, IdealTrappingPotential, axis);
|
||||
if isinf(v), v = NaN; end
|
||||
if isinf(dv), dv = NaN; end
|
||||
IdealTrapFrequency = [v, dv];
|
||||
|
||||
% Extract trap frequencies from the actual trapping potential
|
||||
if options.extract_trap_frequencies
|
||||
[v, dv, popt, pcov] = extractTrapFrequency(Positions, TrappingPotential, axis);
|
||||
if isinf(v), v = NaN; end
|
||||
if isinf(dv), dv = NaN; end
|
||||
TrapFrequency = [v, dv];
|
||||
ExtractedTrapFrequencies = {IdealTrapFrequency, TrapFrequency};
|
||||
else
|
||||
ExtractedTrapFrequencies = {IdealTrapFrequency};
|
||||
end
|
||||
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user