From 8327652883fe6eae8cb6636561f8c074a130b3da Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sun, 11 Jul 2021 14:46:19 +0200 Subject: [PATCH] Shifted some calculations in to other functions. --- .../@MOTSimulator/runSimulation.m | 46 ++++--------------- 1 file changed, 8 insertions(+), 38 deletions(-) diff --git a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m b/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m index 49135f7..e3fdad4 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m +++ b/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m @@ -1,36 +1,6 @@ -function [LoadingRate, StandardError] = runSimulation(this) - n = this.NumberOfAtoms; - %% Calculate Background Collision Time --> Calculate Capture velocity --> Introduce velocity cutoff --> Calculate capture fraction +function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this) - this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); % Make 5% larger than the numerically obtained CV - 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 for initial positions and velocities % - sampling the position distribution this.InitialPositions = this.initialPositionSampling(); % - sampling the velocity distribution @@ -38,18 +8,18 @@ function [LoadingRate, StandardError] = runSimulation(this) %% Solve ODE 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 - FinalDynamicalQuantities = zeros(n,9); + ParticleDynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6); Positions = this.InitialPositions; Velocities = this.InitialVelocities; - parfor Index = 1:n - ret = this.solver(Positions(Index,:), Velocities(Index,:)); - FinalDynamicalQuantities(Index,:) = ret(end,:); + parfor Index = 1:this.NumberOfAtoms + ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:)); progressbar.PB_iterate(); end clear Index + %% Calculate the Loading Rate - [LoadingRate, StandardError] = this.calculateLoadingRate(FinalDynamicalQuantities); + [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ParticleDynamicalQuantities); end \ No newline at end of file