%% 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