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

Major changes - rewrote how autocorrelation, resampling and determination of loading rate is done.

Browse Source
This commit is contained in:
Karthik 2021-07-11 14:47:13 +02:00
parent 8327652883
commit 5b91c20246
1 changed files with 33 additions and 79 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

110
MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m
View File

@ -1,84 +1,38 @@
function [LoadingRate, StandardError] = calculateLoadingRate(this, FinalDynamicalQuantities)
function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities)
switch this.SimulationMode
case "2D"
n = this.NumberOfAtoms;
NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep);
NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps);
TimeCounts = zeros(1, n);
CollisionEvents = zeros(1, n);
NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms);
AutocorrelationFunction = zeros(1, this.NumberOfAtoms);
for i = 1:this.NumberOfAtoms
FinalPosition = FinalDynamicalQuantities(i,1:3);
DivergenceAngle = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2)));
if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0)
if i == 1
NumberOfLoadedAtoms(1) = 1;
else
NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1) + 1;
% Include the stochastic process of background collisions
for AtomIndex = 1:n
TimeCounts(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
end
this.TimeSpentInInteractionRegion = mean(TimeCounts);
for AtomIndex = 1:n
CollisionEvents(AtomIndex) = this.computeCollisionProbability();
end
else
if i == 1
NumberOfLoadedAtoms(1) = 0;
else
NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1);
% 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(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
Velocity = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
if this.exitCondition(Position, Velocity, CollisionEvents(AtomIndex))
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
end
end
end
end
[LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(NumberOfLoadedAtoms);
case "3D"
% Development In progress
end
for i = 1:this.NumberOfAtoms-1
MeanTrappingEfficiency = 0;
MeanLoadingRateShifted = 0;
for j = 1:this.NumberOfAtoms-i
MeanTrappingEfficiency = MeanTrappingEfficiency + NumberOfLoadedAtoms(j)/j;
MeanLoadingRateShifted = MeanLoadingRateShifted + (NumberOfLoadedAtoms(i+j))/(i+j);
AutocorrelationFunction(i) = AutocorrelationFunction(i) + ((NumberOfLoadedAtoms(j)/j).*(NumberOfLoadedAtoms(i+j)/(i+j)));
end
AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanTrappingEfficiency * MeanLoadingRateShifted));
end
if AutocorrelationFunction(1)~=0 % To handle cases where there is no atom loading
AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1);
x = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms);
[FitObject,~] = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100);
CorrelationFactor = FitObject.n; % n is the autocorrelation factor
SumOfAllMeanTrappingEfficiencies = 0;
NumberOfBlocks = floor(this.NumberOfAtoms/(2*CorrelationFactor+1));
MeanTrappingEfficiencyInEachBlock = zeros(1,NumberOfBlocks);
BlockNumberLimit = min(NumberOfBlocks-1,5);
% Jackknife method is a systematic way of obtaining the “standard deviation” error of a
% set of stochastic measurements:
% 1. Calculate average (or some function f) from the full dataset.
% 2. Divide data into M blocks, with block length ≫ correlation factor . This is done in order to
% get rid of autocorrelations; if there are no correlations, block length can be 1.
% 3. For each m = 1 . . .M, take away block m and calculate the average using
% the data from all other blocks.
% 4. Estimate the error by calculating the deviation of the average in each block from average of the full dataset
%Jackknife estimate of the loading rate
for i = 1:NumberOfBlocks-BlockNumberLimit
MeanTrappingEfficiencyInEachBlock(i) = NumberOfLoadedAtoms(this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1))) / (this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1)));
SumOfAllMeanTrappingEfficiencies = SumOfAllMeanTrappingEfficiencies + MeanTrappingEfficiencyInEachBlock(i);
end
MeanTrappingEfficiency = SumOfAllMeanTrappingEfficiencies / (NumberOfBlocks-BlockNumberLimit);
c = integral(@(velocity) sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ...
* velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* this.OvenTemperatureinKelvin))), 0, this.VelocityCutoff);
LoadingRate = MeanTrappingEfficiency * c * this.calculateFreeMolecularRegimeFlux();
% Jackknife estimate of the variance of the loading rate
Variance = 0;
for i = 1:NumberOfBlocks-BlockNumberLimit
Variance = Variance + (MeanTrappingEfficiencyInEachBlock(i) - MeanTrappingEfficiency)^2;
end
StandardError = sqrt((((NumberOfBlocks-BlockNumberLimit) - 1) / (NumberOfBlocks-BlockNumberLimit)) * Variance);
else
LoadingRate = 0;
StandardError = 0;
end
end