diff --git a/Estimations/EstimatesForHybridTrap.m b/Estimations/EstimatesForHybridTrap.m deleted file mode 100644 index 7893eb9..0000000 --- a/Estimations/EstimatesForHybridTrap.m +++ /dev/null @@ -1,368 +0,0 @@ -%% 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