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

More cosmetic changes, corrections to some calculations like the one of flux, re-wrote how one and two parameter scans are done with functions, added a save option to save the entire simulator object with its properties and their values along with the results.

Browse Source
This commit is contained in:
Karthik 2021-07-04 20:24:38 +02:00
parent add5d02ae3
commit d0b4d0b9b8
14 changed files with 293 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m Normal file
View File

@ -0,0 +1,62 @@
function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
p = inputParser;
p.KeepUnmatched = true;
addRequired(p, 'ParameterArray', @isvector)
addRequired(p, 'QuantityOfInterestArray', @ismatrix)
addParameter(p, 'RescalingFactorForParameter', 1, @isscalar)
addParameter(p, 'XLabelString', 'X parameter', @ischar)
addParameter(p, 'ErrorsForYQuantity', false, @islogical)
addParameter(p, 'ErrorsArray', [], @isvector)
addParameter(p, 'RescalingFactorForYQuantity', 1, @isscalar)
addParameter(p, 'YLabelString', 'Y parameter', @ischar)
addParameter(p, 'TitleString', 'One-Parameter Scan', @ischar)
p.parse(XParameter, YQuantity, varargin{:})
XParameter = p.Results.ParameterArray;
RescalingFactorForXParameter = p.Results.RescalingFactorForParameter;
XLabelString = p.Results.XLabelString;
YQuantity = p.Results.QuantityOfInterestArray;
ErrorsForYQuantity = p.Results.ErrorsForYQuantity;
ErrorsArray = p.Results.ErrorsArray;
RescalingFactorForYQuantity = p.Results.RescalingFactorForYQuantity;
YLabelString = p.Results.YLabelString;
TitleString = p.Results.TitleString;
f_h = Helper.getFigureByTag('One-Parameter Scan');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'One-Parameter Scan';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
RescaledXParameter = XParameter .* RescalingFactorForXParameter;
RescaledYQuantity = YQuantity .* RescalingFactorForYQuantity;
if ErrorsForYQuantity
RescaledErrorsArray = ErrorsArray .* RescalingFactorForYQuantity;
errorbar(RescaledXParameter, RescaledYQuantity, RescaledErrorsArray)
else
plot(RescaledXParameter, RescaledYQuantity, '--o', 'Linewidth', 1.5);
end
hXLabel = xlabel(XLabelString);
hYLabel = ylabel(YLabelString);
hTitle = sgtitle(TitleString);
set([hXLabel, hYLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

7
MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m
View File

@ -11,6 +11,7 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg
addParameter(p, 'XLabelString', 'X parameter', @ischar) addParameter(p, 'XLabelString', 'X parameter', @ischar)
addParameter(p, 'RescalingFactorForSecondParameter', 1, @isscalar) addParameter(p, 'RescalingFactorForSecondParameter', 1, @isscalar)
addParameter(p, 'YLabelString', 'Y parameter', @ischar) addParameter(p, 'YLabelString', 'Y parameter', @ischar)
addParameter(p, 'RescalingFactorForQuantityOfInterest', 1, @isscalar)
addParameter(p, 'ZLabelString', 'Z parameter', @ischar) addParameter(p, 'ZLabelString', 'Z parameter', @ischar)
addParameter(p, 'TitleString', 'Two-Parameter Scan', @ischar) addParameter(p, 'TitleString', 'Two-Parameter Scan', @ischar)
@ -23,6 +24,7 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg
RescalingFactorForYParameter = p.Results.RescalingFactorForSecondParameter; RescalingFactorForYParameter = p.Results.RescalingFactorForSecondParameter;
YLabelString = p.Results.YLabelString; YLabelString = p.Results.YLabelString;
ZQuantity = p.Results.QuantityOfInterestArray; ZQuantity = p.Results.QuantityOfInterestArray;
RescalingFactorForZQuantity = p.Results.RescalingFactorForQuantityOfInterest;
ZLabelString = p.Results.ZLabelString; ZLabelString = p.Results.ZLabelString;
TitleString = p.Results.TitleString; TitleString = p.Results.TitleString;
@ -41,12 +43,13 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg
RescaledXParameter = XParameter .* RescalingFactorForXParameter; RescaledXParameter = XParameter .* RescalingFactorForXParameter;
RescaledYParameter = YParameter .* RescalingFactorForYParameter; RescaledYParameter = YParameter .* RescalingFactorForYParameter;
RescaledZQuantity = ZQuantity .* RescalingFactorForZQuantity;
imagesc(RescaledXParameter, RescaledYParameter, ZQuantity(:,:)'); imagesc(RescaledXParameter, RescaledYParameter, RescaledZQuantity(:,:)');
set(gca,'YDir','normal'); set(gca,'YDir','normal');
caxis([min(min(min(ZQuantity))) max(max(max(ZQuantity)))]); caxis([min(min(min(RescaledZQuantity))) max(max(max(RescaledZQuantity)))]);
hXLabel = xlabel(XLabelString); hXLabel = xlabel(XLabelString);
hYLabel = ylabel(YLabelString); hYLabel = ylabel(YLabelString);

28
MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m
View File

@ -70,7 +70,6 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
CaptureVelocity; CaptureVelocity;
VelocityCutoff; VelocityCutoff;
CaptureFraction;
ClausingFactor; ClausingFactor;
ThetaArray; ThetaArray;
AngularDistribution; AngularDistribution;
@ -85,6 +84,9 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
DebugMode; DebugMode;
DoSave; DoSave;
SaveDirectory;
Results;
end end
properties (SetAccess = private, GetAccess = public) properties (SetAccess = private, GetAccess = public)
@ -128,6 +130,8 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
@islogical); @islogical);
addParameter(p, 'SaveData', false,... addParameter(p, 'SaveData', false,...
@islogical); @islogical);
addParameter(p, 'SaveDirectory', pwd,...
@ischar);
p.parse(varargin{:}); p.parse(varargin{:});
@ -144,7 +148,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
s.DebugMode = p.Results.DebugMode; s.DebugMode = p.Results.DebugMode;
s.DoSave = p.Results.SaveData; s.DoSave = p.Results.SaveData;
s.SaveDirectory = p.Results.SaveDirectory;
s.reinitializeSimulator(); s.reinitializeSimulator();
@ -519,12 +523,6 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.VelocityCutoff(this) function ret = get.VelocityCutoff(this)
ret = this.VelocityCutoff; ret = this.VelocityCutoff;
end end
function set.CaptureFraction(this,val)
this.CaptureFraction = val;
end
function ret = get.CaptureFraction(this)
ret = this.CaptureFraction;
end
function set.ClausingFactor(this,val) function set.ClausingFactor(this,val)
this.ClausingFactor = val; this.ClausingFactor = val;
end end
@ -570,6 +568,20 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.DoSave(this) function ret = get.DoSave(this)
ret = this.DoSave; ret = this.DoSave;
end end
function set.SaveDirectory(this, val)
this.SaveDirectory = val;
end
function ret = get.SaveDirectory(this)
ret = this.SaveDirectory;
end
function set.Results(this, val)
this.Results = val;
end
function ret = get.Results(this)
ret = this.Results;
end
end % - setters and getters end % - setters and getters
methods methods

2
MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m
View File

@ -11,7 +11,7 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
Detuning = this.PushBeamDetuning - DopplerShift; Detuning = this.PushBeamDetuning - DopplerShift;
a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveVector * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* ... a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveVector * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* ...
(SaturationIntensity/(1 + SaturationIntensity + 4 * (Detuning./this.PushBeamLinewidth).^2)); (SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2)));
if this.SpontaneousEmission if this.SpontaneousEmission
a_scatter = this.accelerationDueToSpontaneousEmissionProcess(SaturationIntensity, SaturationIntensity, Detuning, this.PushBeamLinewidth, this.PushBeamWaveVector); a_scatter = this.accelerationDueToSpontaneousEmissionProcess(SaturationIntensity, SaturationIntensity, Detuning, this.PushBeamLinewidth, this.PushBeamWaveVector);

2
MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m
View File

@ -2,7 +2,7 @@ function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationInten
Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1]; Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1];
Vector = Vector./norm(Vector); Vector = Vector./norm(Vector);
ScatteringRate = 0.5 * SaturationIntensity * Linewidth / (1 + TotalSaturationIntensity + 4 * (Detuning/Linewidth)^2); ScatteringRate = (0.5 * SaturationIntensity * Linewidth) / (1 + TotalSaturationIntensity + (4 * (Detuning/Linewidth)^2));
NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep); NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep);
if NumberOfScatteringEvents > 0 if NumberOfScatteringEvents > 0

1
MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m
View File

@ -20,6 +20,5 @@ function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
break; break;
end end
end end
clear Index
end end

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

@ -1,4 +1,4 @@
function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFraction, ClausingFactor, FinalDynamicalQuantities) function [LoadingRate, StandardError] = calculateLoadingRate(this, FinalDynamicalQuantities)
NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms); NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms);
AutocorrelationFunction = zeros(1, this.NumberOfAtoms); AutocorrelationFunction = zeros(1, this.NumberOfAtoms);
@ -23,41 +23,58 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti
end end
for i = 1:this.NumberOfAtoms-1 for i = 1:this.NumberOfAtoms-1
MeanLoadingRate = 0; MeanTrappingEfficiency = 0;
MeanLoadingRateShifted = 0; MeanLoadingRateShifted = 0;
for j = 1:this.NumberOfAtoms-i for j = 1:this.NumberOfAtoms-i
MeanLoadingRate = MeanLoadingRate + NumberOfLoadedAtoms(j)/j; MeanTrappingEfficiency = MeanTrappingEfficiency + NumberOfLoadedAtoms(j)/j;
MeanLoadingRateShifted = MeanLoadingRateShifted + (NumberOfLoadedAtoms(i+j))/(i+j); MeanLoadingRateShifted = MeanLoadingRateShifted + (NumberOfLoadedAtoms(i+j))/(i+j);
AutocorrelationFunction(i) = AutocorrelationFunction(i) + ((NumberOfLoadedAtoms(j)/j).*(NumberOfLoadedAtoms(i+j)/(i+j))); AutocorrelationFunction(i) = AutocorrelationFunction(i) + ((NumberOfLoadedAtoms(j)/j).*(NumberOfLoadedAtoms(i+j)/(i+j)));
end end
AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanLoadingRate * MeanLoadingRateShifted)); AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanTrappingEfficiency * MeanLoadingRateShifted));
end end
if AutocorrelationFunction(1)~=0 % In case no atom loading if AutocorrelationFunction(1)~=0 % To handle cases where there is no atom loading
AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1); AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1);
x = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms); x = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms);
[FitObject,~] = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100); [FitObject,~] = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100);
n = FitObject.n; % n is the autocorrelation factor CorrelationFactor = FitObject.n; % n is the autocorrelation factor
MeanLoadingRate = 0; SumOfAllMeanTrappingEfficiencies = 0;
NumberOfBins = floor(this.NumberOfAtoms/(2*n+1)); NumberOfBlocks = floor(this.NumberOfAtoms/(2*CorrelationFactor+1));
LoadingRateError = zeros(1,NumberOfBins); MeanTrappingEfficiencyInEachBlock = zeros(1,NumberOfBlocks);
BinNumberLimit = min(NumberOfBins-1,5); BlockNumberLimit = min(NumberOfBlocks-1,5);
for i = 1:NumberOfBins-BinNumberLimit
LoadingRateError(i) = NumberOfLoadedAtoms(this.NumberOfAtoms-ceil((i-1)*(2*n+1))) / ... % Jackknife method is a systematic way of obtaining the “standard deviation” error of a
(this.NumberOfAtoms-ceil((i-1)*(2*n+1))); % set of stochastic measurements:
MeanLoadingRate = MeanLoadingRate + LoadingRateError(i); % 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 end
MeanLoadingRate = MeanLoadingRate /(NumberOfBins-BinNumberLimit); MeanTrappingEfficiency = SumOfAllMeanTrappingEfficiencies / (NumberOfBlocks-BlockNumberLimit);
StandardError = 0; c = integral(@(velocity) sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ...
for i = 1:NumberOfBins-BinNumberLimit * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
StandardError = StandardError + (MeanLoadingRate-LoadingRateError(i))^2; * 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 end
StandardError = sqrt(StandardError/(NumberOfBins-BinNumberLimit));
LoadingRate = (MeanLoadingRate * this.calculateFreeMolecularRegimeFlux() * CaptureFraction * ClausingFactor); StandardError = sqrt((((NumberOfBlocks-BlockNumberLimit) - 1) / (NumberOfBlocks-BlockNumberLimit)) * Variance);
else else
LoadingRate = 0; LoadingRate = 0;

4
MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m
View File

@ -6,9 +6,7 @@ function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVecto
%One side of parallelogram = PositionVectorFromWaveVectorOrigin %One side of parallelogram = PositionVectorFromWaveVectorOrigin
%Base = Wavevector %Base = Wavevector
%Area = One side of parallelogram X Base %Area = One side of parallelogram X Base
DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector); % Slow DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector);
%DistanceBetweenAtomAndLaserBeamAxis = norm((WaveVector*WaveVector')*PositionVectorFromWaveVectorOrigin-(WaveVector*PositionVectorFromWaveVectorOrigin')*WaveVector)./ ...
% (WaveVector(1)^2+WaveVector(2)^2+WaveVector(3)^2); % Faster
if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius
ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2); ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2);

63
MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m Normal file
View File

@ -0,0 +1,63 @@
function [LoadingRateArray, StandardErrorArray] = doOneParameterScan(this, ParameterName, ParameterArray, varargin)
p = inputParser;
p.KeepUnmatched = true;
addRequired(p, 'ClassObject' , @isobject);
addRequired(p, 'ParameterName' , @ischar);
addRequired(p, 'ParameterArray' , @isvector);
addParameter(p, 'ChangeRelatedParameter', false, @islogical);
addParameter(p, 'RelatedParameterName', 'none', @ischar);
addParameter(p, 'RelatedParameterArray', length(ParameterArray), @isvector);
addParameter(p, 'ChangeInitialConditions', false, @islogical);
addParameter(p, 'ParameterNameArray', {}, @iscell);
addParameter(p, 'ParameterValueArray', {}, @iscell);
p.parse(this, ParameterName, ParameterArray, varargin{:})
ParameterName = p.Results.ParameterName;
ParameterArray = p.Results.ParameterArray;
ChangeRelatedParameter = p.Results.ChangeRelatedParameter;
RelatedParameterName = p.Results.RelatedParameterName;
RelatedParameterArray = p.Results.RelatedParameterArray;
ChangeInitialConditions = p.Results.ChangeInitialConditions;
ParameterNameArray = p.Results.ParameterNameArray;
ParameterValueArray = p.Results.ParameterValueArray;
NumberOfPointsForParam = length(ParameterArray);
LoadingRateArray = zeros(1,NumberOfPointsForParam);
StandardErrorArray = zeros(1,NumberOfPointsForParam);
for i=1:NumberOfPointsForParam
eval(sprintf('OptionsStruct.%s = %d;', ParameterName, ParameterArray(i)));
if ChangeRelatedParameter
eval(sprintf('OptionsStruct.%s = %d;', RelatedParameterName, RelatedParameterArray(i)));
end
if ChangeInitialConditions
for j = 1:length(ParameterNameArray)
eval(sprintf('OptionsStruct.%s = %d;', ParameterNameArray{j}, ParameterValueArray{j}));
end
end
options = Helper.convertstruct2cell(OptionsStruct);
this.setInitialConditions(options{:});
tic
[LoadingRate, StandardError] = this.runSimulation();
LoadingRateArray(i) = LoadingRate;
StandardErrorArray(i) = StandardError;
end
if this.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
LoadingRate.Errors = StandardErrorArray;
this.Results = LoadingRate;
SaveFolder = [this.SaveDirectory filesep 'Results'];
Filename = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = this;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
end

30
MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m
View File

@ -1,4 +1,4 @@
function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ... function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ...
SecondParameterName, SecondParameterArray, varargin) SecondParameterName, SecondParameterArray, varargin)
p = inputParser; p = inputParser;
@ -15,6 +15,9 @@ function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstPa
addParameter(p, 'RelatedParameterName', 'none', @ischar); addParameter(p, 'RelatedParameterName', 'none', @ischar);
addParameter(p, 'RelatedParameterArray', length(FirstParameterArray), @isvector); addParameter(p, 'RelatedParameterArray', length(FirstParameterArray), @isvector);
addParameter(p, 'ChangeInitialConditions', false, @islogical);
addParameter(p, 'ParameterNameArray', {}, @iscell);
addParameter(p, 'ParameterValueArray', {}, @iscell);
p.parse(this, FirstParameterName, FirstParameterArray, ... p.parse(this, FirstParameterName, FirstParameterArray, ...
SecondParameterName, SecondParameterArray, varargin{:}) SecondParameterName, SecondParameterArray, varargin{:})
@ -27,9 +30,14 @@ function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstPa
Order = p.Results.Order; Order = p.Results.Order;
RelatedParameterName = p.Results.RelatedParameterName; RelatedParameterName = p.Results.RelatedParameterName;
RelatedParameterArray = p.Results.RelatedParameterArray; RelatedParameterArray = p.Results.RelatedParameterArray;
ChangeInitialConditions = p.Results.ChangeInitialConditions;
ParameterNameArray = p.Results.ParameterNameArray;
ParameterValueArray = p.Results.ParameterValueArray;
NumberOfPointsForFirstParam = length(FirstParameterArray); NumberOfPointsForFirstParam = length(FirstParameterArray);
NumberOfPointsForSecondParam = length(SecondParameterArray); NumberOfPointsForSecondParam = length(SecondParameterArray);
LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
for i=1:NumberOfPointsForFirstParam for i=1:NumberOfPointsForFirstParam
eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i))); eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i)));
@ -42,11 +50,29 @@ function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstPa
if ChangeRelatedParameter && Order == 2 if ChangeRelatedParameter && Order == 2
eval(sprintf('OptionsStruct.%s = %d;', RelatedParameterName, RelatedParameterArray(j))); eval(sprintf('OptionsStruct.%s = %d;', RelatedParameterName, RelatedParameterArray(j)));
end end
if ChangeInitialConditions
for k = 1:length(ParameterNameArray)
eval(sprintf('OptionsStruct.%s = %d;', ParameterNameArray{k}, ParameterValueArray{k}));
end
end
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
this.setInitialConditions(options{:}); this.setInitialConditions(options{:});
tic tic
[LoadingRate, ~] = this.runSimulation(); [LoadingRate, StandardError] = this.runSimulation();
LoadingRateArray(i,j) = LoadingRate; LoadingRateArray(i,j) = LoadingRate;
StandardErrorArray(i,j) = StandardError;
end end
end end
if this.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
LoadingRate.Errors = StandardErrorArray;
this.Results = LoadingRate;
SaveFolder = [this.SaveDirectory filesep 'Results'];
Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = this;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
end end

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

@ -5,13 +5,6 @@ function [LoadingRate, StandardError] = runSimulation(this)
this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); % Make 5% larger than the numerically obtained CV 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 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. %velocity is only along the x-axis, we take the first term which is the x-component of the velocity.
VelocityDistributionFunction = @(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)));
this.CaptureFraction = integral(VelocityDistributionFunction, 0, this.VelocityCutoff) / integral(VelocityDistributionFunction, 0, Inf);
%% Calculate the Clausing Factor %% Calculate the Clausing Factor
% Compute the angular distribution of the atomic beam % Compute the angular distribution of the atomic beam
@ -58,5 +51,5 @@ function [LoadingRate, StandardError] = runSimulation(this)
end end
clear Index clear Index
%% Calculate the Loading Rate %% Calculate the Loading Rate
[LoadingRate, StandardError] = this.calculateLoadingRate(this.CaptureFraction, this.ClausingFactor, FinalDynamicalQuantities); [LoadingRate, StandardError] = this.calculateLoadingRate(FinalDynamicalQuantities);
end end

1
MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m
View File

@ -66,6 +66,7 @@ function setInitialConditions(this, varargin)
this.CoolingBeamWaveVector = this.BlueWaveVector; this.CoolingBeamWaveVector = this.BlueWaveVector;
this.CoolingBeamDetuning = this.BlueDetuning; this.CoolingBeamDetuning = this.BlueDetuning;
this.CoolingBeamRadius = this.BlueBeamRadius; this.CoolingBeamRadius = this.BlueBeamRadius;
this.CoolingBeamWaist = this.BlueBeamWaist;
this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity; this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity;
this.SidebandBeamRadius = this.BlueBeamRadius; this.SidebandBeamRadius = this.BlueBeamRadius;
this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity; this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity;

89
MOT Capture Process Simulation/test_MOTSimulator.m
View File

@ -6,13 +6,14 @@ OptionsStruct = struct;
OptionsStruct.SimulationMode = '2D'; OptionsStruct.SimulationMode = '2D';
OptionsStruct.TimeStep = 50e-06; % in s OptionsStruct.TimeStep = 50e-06; % in s
OptionsStruct.SimulationTime = 4e-03; % in s OptionsStruct.SimulationTime = 4e-03; % in s
OptionsStruct.SpontaneousEmission = false; OptionsStruct.SpontaneousEmission = true;
OptionsStruct.Sideband = false; OptionsStruct.Sideband = true;
OptionsStruct.ZeemanSlowerBeam = false; OptionsStruct.ZeemanSlowerBeam = false;
OptionsStruct.Gravity = true; OptionsStruct.Gravity = true;
OptionsStruct.AtomicBeamCollision = true; OptionsStruct.AtomicBeamCollision = true;
OptionsStruct.DebugMode = false; OptionsStruct.DebugMode = false;
OptionsStruct.SaveData = false; OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation';
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
Simulator = MOTSimulator(options{:}); Simulator = MOTSimulator(options{:});
@ -52,7 +53,7 @@ Plotting.visualizeMagneticField(Simulator, XAxisRange, YAxisRange, ZAxisRange)
%% - Plot MFP & VP for different temperatures %% - Plot MFP & VP for different temperatures
TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius
Plotting.plotMeanFreePathAndVapourPressure(TemperatureinCelsius) Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
%% - Plot the Free Molecular Flux for different temperatures %% - Plot the Free Molecular Flux for different temperatures
Temperature = [900, 950, 1000, 1050, 1100]; % Temperature' Temperature = [900, 950, 1000, 1050, 1100]; % Temperature'
@ -82,45 +83,54 @@ IncidentAtomPosition = 0;
Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition) Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
%% - Scan parameters %% - Scan parameters
%Use a for loop to change the parameter during set initial conditions
%Run simulation
% TWO-PARAMETER SCAN % ONE-PARAMETER SCAN
NumberOfPointsForFirstParam = 2; %iterations of the simulation NumberOfPointsForParam = 20; %iterations of the simulation
NumberOfPointsForSecondParam = 2; % Scan Cooling Beam Power
LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); PowerArray = linspace(0.1, 1.0, NumberOfPointsForParam) * Simulator.TotalPower;
% Scan Cooling Beam Detuning
% DetuningArray = linspace(-0.5,-10, NumberOfPointsForParam) * Helper.PhysicsConstants.BlueLinewidth;
% Scan Sideband Detuning and Power Ratio OptionsStruct = struct;
DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam); OptionsStruct.ChangeInitialConditions = true;
PowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam); OptionsStruct.ParameterNameArray = {'NumberOfAtoms'};
OptionsStruct.ParameterValueArray = {10000};
OptionsStruct.NumberOfAtoms = 5000; options = Helper.convertstruct2cell(OptionsStruct);
tStart = tic; tStart = tic;
for i=1:NumberOfPointsForFirstParam [LoadingRateArray, ~] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
OptionsStruct.SidebandDetuning = DetuningArray(i) * Helper.PhysicsConstants.BlueLinewidth;
for j=1:NumberOfPointsForSecondParam
OptionsStruct.BluePower = PowerArray(j) * Simulator.TotalPower;
OptionsStruct.SidebandPower = Simulator.TotalPower - OptionsStruct.BluePower;
options = Helper.convertstruct2cell(OptionsStruct);
Simulator.setInitialConditions(options{:});
tic
[LoadingRate, ~] = Simulator.runSimulation();
LoadingRateArray(i,j) = LoadingRate;
end
end
tEnd = toc(tStart); tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
%% TWO-PARAMETER SCAN FUNCTION clear OptionsStruct
NumberOfPointsForFirstParam = 2; %iterations of the simulation % - Plot results
NumberOfPointsForSecondParam = 2;
Simulator.NumberOfAtoms = 5000; ParameterArray = PowerArray .* Simulator.TotalPower;
QuantityOfInterestArray = LoadingRateArray;
OptionsStruct = struct;
OptionsStruct.RescalingFactorForParameter = 1000;
OptionsStruct.XLabelString = 'Cooling Beam Power (mW)';
OptionsStruct.RescalingFactorForYQuantity = 1e-16;
OptionsStruct.ErrorsForYQuantity = false;
OptionsStruct.YLabelString = 'Loading rate (x 10^{16} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotting.plotResultForOneParameterScan(ParameterArray, QuantityOfInterestArray, options{:})
clear OptionsStruct
%% TWO-PARAMETER SCAN
NumberOfPointsForParam = 20; %iterations of the simulation
NumberOfPointsForSecondParam = 10;
% Scan Sideband Detuning and Power Ratio % Scan Sideband Detuning and Power Ratio
DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; DetuningArray = linspace(-0.5,-6, NumberOfPointsForParam) * Helper.PhysicsConstants.BlueLinewidth;
SidebandPowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam) * Simulator.TotalPower; SidebandPowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam) * Simulator.TotalPower;
BluePowerArray = Simulator.TotalPower - SidebandPowerArray; BluePowerArray = Simulator.TotalPower - SidebandPowerArray;
@ -129,18 +139,22 @@ OptionsStruct.ChangeRelatedParameter = true;
OptionsStruct.Order = 2; %Change after first parameter = 1, Change after second parameter = 2 OptionsStruct.Order = 2; %Change after first parameter = 1, Change after second parameter = 2
OptionsStruct.RelatedParameterName = 'BluePower'; OptionsStruct.RelatedParameterName = 'BluePower';
OptionsStruct.RelatedParameterArray = BluePowerArray; OptionsStruct.RelatedParameterArray = BluePowerArray;
OptionsStruct.ChangeInitialConditions = true;
OptionsStruct.ParameterNameArray = {'NumberOfAtoms'};
OptionsStruct.ParameterValueArray = {10000};
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
tStart = tic; tStart = tic;
LoadingRateArray = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:}); [LoadingRateArray, StandardErrorArray] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
tEnd = toc(tStart); tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
clear OptionsStruct clear OptionsStruct
%% - Plot results
FirstParameterArray = DetuningArray * Helper.PhysicsConstants.BlueLinewidth; % - Plot results
SecondParameterArray = (Simulator.TotalPower - (PowerArray * Simulator.TotalPower));
FirstParameterArray = DetuningArray;
SecondParameterArray = SidebandPowerArray;
QuantityOfInterestArray = LoadingRateArray; QuantityOfInterestArray = LoadingRateArray;
OptionsStruct = struct; OptionsStruct = struct;
@ -148,7 +162,8 @@ OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.B
OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)'; OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000; OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Sideband Power (mW)'; OptionsStruct.YLabelString = 'Sideband Power (mW)';
OptionsStruct.ZLabelString = 'Loading rate (atoms/s)'; OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-16;
OptionsStruct.ZLabelString = 'Loading rate (x 10^{16} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100); OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);