66 lines
2.9 KiB
Mathematica
66 lines
2.9 KiB
Mathematica
|
function [v, dv] = extractTrapFrequency(Positions, TrappingPotential, options)
|
||
|
|
% 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.
|
||
|
|
|
||
|
|
axis = options.Axis;
|
||
|
|
|
||
|
|
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);
|
||
|
|
|
||
|
|
% Perform curve fitting using harmonic potential
|
||
|
|
try
|
||
|
|
% Define the fit model, treating 'm' as a fixed constant
|
||
|
|
harmonic_potential = fittype(@(v, xoffset, yoffset, m, pos) Potentials.generateHarmonicPotential(pos, v, xoffset, yoffset, m), ...
|
||
|
|
'independent', 'pos', ...
|
||
|
|
'coefficients', {'v', 'xoffset', 'yoffset'}, ...
|
||
|
|
'problem', 'm');
|
||
|
|
|
||
|
|
% Fit the potential data to the harmonic potential model
|
||
|
|
fitObject = fit(xdata(:), Potential(:), harmonic_potential, 'StartPoint', [1e3, tmp_pos(center_idx), -100], 'problem', options.Mass);
|
||
|
|
|
||
|
|
% Get 95% confidence intervals (default)
|
||
|
|
ci = confint(fitObject);
|
||
|
|
|
||
|
|
% Extract the parameter values and their uncertainties
|
||
|
|
params = coeffvalues(fitObject); % Extract the fitted parameter values
|
||
|
|
lower_bound = ci(1, :); % Lower bound of the confidence intervals
|
||
|
|
upper_bound = ci(2, :); % Upper bound of the confidence intervals
|
||
|
|
|
||
|
|
% Compute the uncertainties (standard error can be estimated as half the width of the CI)
|
||
|
|
uncertainties = (upper_bound - lower_bound) / 2;
|
||
|
|
|
||
|
|
v = fitObject.f; % Extract fitted frequency
|
||
|
|
dv = NaN; % Uncertainty in the frequency
|
||
|
|
catch
|
||
|
|
% In case of fitting failure, return NaNs
|
||
|
|
v = NaN;
|
||
|
|
dv = NaN;
|
||
|
|
end
|
||
|
|
end
|