From c19a5140277cc9e36aa874d1e8e1b139ea5a2f73 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 16 Jul 2021 15:57:16 +0200 Subject: [PATCH] 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. --- .../@TwoDimensionalMOT/calculateLoadingRate.m | 51 ++++++++++++++----- .../jackknifeErrorEstimation.m | 50 ++++++++++++++++++ 2 files changed, 88 insertions(+), 13 deletions(-) create mode 100644 +Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m diff --git a/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m b/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m index ee42339..eada4cb 100644 --- a/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m +++ b/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m @@ -1,9 +1,9 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ovenObj) n = this.NumberOfAtoms; DynamicalQuantities = this.ParticleDynamicalQuantities; - NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep); - NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps); CollisionEvents = zeros(1, n); + NumberOfLoadedAtoms = zeros(1, n); + % Include the stochastic process of background collisions for AtomIndex = 1:n @@ -12,17 +12,42 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate end % 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; + + 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 - 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 - - [LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(ovenObj, NumberOfLoadedAtoms); end diff --git a/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m b/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m new file mode 100644 index 0000000..abca489 --- /dev/null +++ b/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m @@ -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 \ No newline at end of file