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

Shifted some calculations in to other functions.

Browse Source
This commit is contained in:
Karthik 2021-07-11 14:46:19 +02:00
parent 5e933b0324
commit 8327652883
1 changed files with 8 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
MOT Capture Process Simulation/@MOTSimulator/runSimulation.m
View File

@ -1,36 +1,6 @@
function [LoadingRate, StandardError] = runSimulation(this) function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this)
n = this.NumberOfAtoms;
%% Calculate Background Collision Time --> Calculate Capture velocity --> Introduce velocity cutoff --> Calculate capture fraction
this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); % Make 5% larger than the numerically obtained CV %% - Sampling for initial positions and velocities
this.VelocityCutoff = this.CaptureVelocity(1); %Should be the magnitude of the 3-D velocity vector but since here the obtained capture
%velocity is only along the x-axis, we take the first term which is the x-component of the velocity.
%% Calculate the Clausing Factor
% Compute the angular distribution of the atomic beam
ThetaArray = linspace(0.0001, pi/2, 1000);
AngularDistribution = zeros(1,length(ThetaArray));
parfor i = 1:length(ThetaArray)
AngularDistribution(i) = this.angularDistributionFunction(ThetaArray(i));
end
% Numerically integrate the angular distribution over the full solid angle
NormalizationConstant = 0;
for j = 1:length(ThetaArray)
if ThetaArray(j) <= this.NozzleExitDivergence
NormalizationConstant = NormalizationConstant + (2 * pi * sin(ThetaArray(j)) * AngularDistribution(j) * (ThetaArray(2)-ThetaArray(1)));
end
end
this.ClausingFactor = (1/pi) * NormalizationConstant; %The complete intergration will give pi * ClausingFactor.
%Therefore, the Clausing Factor needs to be extracted from the
%result of the above integration by dividing out pi
this.ThetaArray = ThetaArray;
this.AngularDistribution = AngularDistribution;
this.NormalizationConstantForAngularDistribution = NormalizationConstant;
%%
% - sampling the position distribution % - sampling the position distribution
this.InitialPositions = this.initialPositionSampling(); this.InitialPositions = this.initialPositionSampling();
% - sampling the velocity distribution % - sampling the velocity distribution
@ -38,18 +8,18 @@ function [LoadingRate, StandardError] = runSimulation(this)
%% Solve ODE %% Solve ODE
progressbar = Helper.parforNotifications(); progressbar = Helper.parforNotifications();
progressbar.PB_start(n,'Message',['Simulating capture process for ' num2str(n,'%.0f') ' atoms:']); progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']);
% calculate the final position of the atoms % calculate the final position of the atoms
FinalDynamicalQuantities = zeros(n,9); ParticleDynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
Positions = this.InitialPositions; Positions = this.InitialPositions;
Velocities = this.InitialVelocities; Velocities = this.InitialVelocities;
parfor Index = 1:n parfor Index = 1:this.NumberOfAtoms
ret = this.solver(Positions(Index,:), Velocities(Index,:)); ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
FinalDynamicalQuantities(Index,:) = ret(end,:);
progressbar.PB_iterate(); progressbar.PB_iterate();
end end
clear Index clear Index
%% Calculate the Loading Rate %% Calculate the Loading Rate
[LoadingRate, StandardError] = this.calculateLoadingRate(FinalDynamicalQuantities); [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ParticleDynamicalQuantities);
end end