karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Added back the captured atom counting and the jackknife error estimation method as developed by Jianshun Gao, retained the bootstrap method and the counting that goes with it now corrected to prevent counting the same atoms repeatedly and also to reduce the atom count if they are lost. There is now an agreement between the two methods as to what the loading rate is with the bootstrap method giving a slightly lower estimate.

Browse Source
This commit is contained in:
Karthik 2021-07-16 15:57:16 +02:00
parent 5407da354b
commit c19a514027
2 changed files with 88 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m
View File

@ -1,9 +1,9 @@
function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ovenObj) function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ovenObj)
n = this.NumberOfAtoms; n = this.NumberOfAtoms;
DynamicalQuantities = this.ParticleDynamicalQuantities; DynamicalQuantities = this.ParticleDynamicalQuantities;
NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep);
NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps);
CollisionEvents = zeros(1, n); CollisionEvents = zeros(1, n);
NumberOfLoadedAtoms = zeros(1, n);
% Include the stochastic process of background collisions % Include the stochastic process of background collisions
for AtomIndex = 1:n for AtomIndex = 1:n
@ -12,17 +12,42 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
end end
% Count the number of loaded atoms subject to conditions % Count the number of loaded atoms subject to conditions
for TimeIndex = 1:NumberOfTimeSteps
if TimeIndex ~= 1
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
end
for AtomIndex = 1:n
Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
if this.exitCondition(Position, CollisionEvents(AtomIndex))
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
end
end
end
[LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(ovenObj, NumberOfLoadedAtoms); switch this.ErrorEstimationMethod
case 'bootstrap'
NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep);
NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps);
LoadedAtomIndices = [];
for TimeIndex = 1:NumberOfTimeSteps
if TimeIndex ~= 1
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
end
for AtomIndex = 1:n
Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
if this.exitCondition(Position, CollisionEvents(AtomIndex))
if ~ismember(AtomIndex, LoadedAtomIndices)
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
LoadedAtomIndices(end+1) = AtomIndex;
end
else
if ismember(AtomIndex, LoadedAtomIndices)
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) - 1;
LoadedAtomIndices(LoadedAtomIndices==AtomIndex) = [];
end
end
end
end
[LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(ovenObj, NumberOfLoadedAtoms);
case 'jackknife'
for AtomIndex = 1:n
if AtomIndex ~= 1
NumberOfLoadedAtoms(AtomIndex) = NumberOfLoadedAtoms(AtomIndex-1);
end
Position = squeeze(DynamicalQuantities(AtomIndex, end, 1:3))';
if this.exitCondition(Position, CollisionEvents(AtomIndex))
NumberOfLoadedAtoms(AtomIndex) = NumberOfLoadedAtoms(AtomIndex) + 1;
end
end
[LoadingRate, StandardError, ConfidenceInterval] = jackknifeErrorEstimation(this, ovenObj, NumberOfLoadedAtoms);
end
end end

50
+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m Normal file
View File

@ -0,0 +1,50 @@
function [LoadingRate, StandardError, ConfidenceInterval] = jackknifeErrorEstimation(this, ovenObj, NumberOfLoadedAtoms)
n = this.NumberOfAtoms;
Autocorrelation = zeros(1, n);
for i = 1:n-1
FirstTerm = 0;
SecondTerm = 0;
for j = 1:n-i
FirstTerm = FirstTerm + NumberOfLoadedAtoms(j) / j;
SecondTerm = SecondTerm + (NumberOfLoadedAtoms(i+j)) / (i+j);
Autocorrelation(i) = Autocorrelation(i) + ((NumberOfLoadedAtoms(j) / j) .*(NumberOfLoadedAtoms(i+j) / (i+j)));
end
Autocorrelation(i) = (1/(n-i)) * (Autocorrelation(i) - ((1/(n-i)) * FirstTerm * SecondTerm));
end
if Autocorrelation(1)~=0
Autocorrelation = Autocorrelation./Autocorrelation(1);
x = linspace(1,n,n);
[FitParams,~] = fit(x',Autocorrelation',"exp(-x/tau)", 'Startpoint', 100);
CorrelationFactor = FitParams.tau;
SampleLength = 2*CorrelationFactor+1;
NumberOfJackknifeSamples = floor(n/SampleLength);
CaptureRatioInEachSample = zeros(1,NumberOfJackknifeSamples);
SampleNumberLimit = min(NumberOfJackknifeSamples-1,5);
for i=1:NumberOfJackknifeSamples-SampleNumberLimit
CaptureRatioInEachSample(i) = NumberOfLoadedAtoms(n-ceil((i-1)*SampleLength))/(n-ceil((i-1)*SampleLength));
end
MeanCaptureRatio = sum(CaptureRatioInEachSample) / (NumberOfJackknifeSamples-SampleNumberLimit);
LoadingRate = MeanCaptureRatio * ovenObj.ReducedFlux;
Variance=0;
for i=1:NumberOfJackknifeSamples-SampleNumberLimit
Variance=Variance+(CaptureRatioInEachSample(i) - MeanCaptureRatio)^2;
end
StandardError = sqrt(Variance/(NumberOfJackknifeSamples-SampleNumberLimit));
ConfidenceInterval = LoadingRate + 1.96*StandardError; % 95% Confidence Intervals
else
LoadingRate = nan;
StandardError = nan;
ConfidenceInterval = [nan nan];
end
end