From d0b4d0b9b84f140b55cb75815a6aa6c983b51a4f Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sun, 4 Jul 2021 20:24:38 +0200 Subject: [PATCH] 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. --- .../+Helper/PhysicsConstants.m | 24 ++--- .../+Plotting/plotResultForOneParameterScan.m | 62 +++++++++++++ .../+Plotting/plotResultForTwoParameterScan.m | 7 +- .../@MOTSimulator/MOTSimulator.m | 28 ++++-- .../@MOTSimulator/accelerationDueToPushBeam.m | 2 +- ...elerationDueToSpontaneousEmissionProcess.m | 2 +- .../@MOTSimulator/calculateCaptureVelocity.m | 1 - .../@MOTSimulator/calculateLoadingRate.m | 67 ++++++++------ .../calculateLocalSaturationIntensity.m | 4 +- .../@MOTSimulator/doOneParameterScan.m | 63 +++++++++++++ .../@MOTSimulator/doTwoParameterScan.m | 38 ++++++-- .../@MOTSimulator/runSimulation.m | 9 +- .../@MOTSimulator/setInitialConditions.m | 1 + .../test_MOTSimulator.m | 89 +++++++++++-------- 14 files changed, 293 insertions(+), 104 deletions(-) create mode 100644 MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m diff --git a/MOT Capture Process Simulation/+Helper/PhysicsConstants.m b/MOT Capture Process Simulation/+Helper/PhysicsConstants.m index 8262415..8f8cc49 100644 --- a/MOT Capture Process Simulation/+Helper/PhysicsConstants.m +++ b/MOT Capture Process Simulation/+Helper/PhysicsConstants.m @@ -23,19 +23,19 @@ classdef PhysicsConstants < handle GravitationalAcceleration = 9.80553; % Dy specific constants - Dy164Mass = 163.929174751*1.66053878283E-27; + Dy164Mass = 163.929174751*1.66053878283E-27; Dy164IsotopicAbundance = 0.2826; - BlueWavelength = 421.291e-9; - BlueLandegFactor = 1.22; - BlueLifetime = 4.94e-9; - BlueLinewidth = 2.02e8; - OrangeWavelength = 626.086e-9; - OrangeLandegFactor = 1.29; - OrangeLifetime = 1.2e-6; - OrangeLinewidth = 8.5e5; - PushBeamLifetime = 1.2e-6; - PushBeamWaveLength = 626.086e-9; - PushBeamLinewidth = 8.5e5; + BlueWavelength = 421.291e-9; + BlueLandegFactor = 1.22; + BlueLifetime = 4.94e-9; + BlueLinewidth = 2.02e8; + OrangeWavelength = 626.086e-9; + OrangeLandegFactor = 1.29; + OrangeLifetime = 1.2e-6; + OrangeLinewidth = 8.5e5; + PushBeamLifetime = 1.2e-6; + PushBeamWaveLength = 626.086e-9; + PushBeamLinewidth = 8.5e5; end methods diff --git a/MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m b/MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m new file mode 100644 index 0000000..74ef9a5 --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m @@ -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 \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m b/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m index 91bafd6..71ff77d 100644 --- a/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m +++ b/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m @@ -11,6 +11,7 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg addParameter(p, 'XLabelString', 'X parameter', @ischar) addParameter(p, 'RescalingFactorForSecondParameter', 1, @isscalar) addParameter(p, 'YLabelString', 'Y parameter', @ischar) + addParameter(p, 'RescalingFactorForQuantityOfInterest', 1, @isscalar) addParameter(p, 'ZLabelString', 'Z parameter', @ischar) addParameter(p, 'TitleString', 'Two-Parameter Scan', @ischar) @@ -23,6 +24,7 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg RescalingFactorForYParameter = p.Results.RescalingFactorForSecondParameter; YLabelString = p.Results.YLabelString; ZQuantity = p.Results.QuantityOfInterestArray; + RescalingFactorForZQuantity = p.Results.RescalingFactorForQuantityOfInterest; ZLabelString = p.Results.ZLabelString; TitleString = p.Results.TitleString; @@ -41,12 +43,13 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg RescaledXParameter = XParameter .* RescalingFactorForXParameter; RescaledYParameter = YParameter .* RescalingFactorForYParameter; + RescaledZQuantity = ZQuantity .* RescalingFactorForZQuantity; - imagesc(RescaledXParameter, RescaledYParameter, ZQuantity(:,:)'); + imagesc(RescaledXParameter, RescaledYParameter, RescaledZQuantity(:,:)'); 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); hYLabel = ylabel(YLabelString); diff --git a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m index 9cb189b..a83f68c 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m +++ b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m @@ -70,7 +70,6 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable CaptureVelocity; VelocityCutoff; - CaptureFraction; ClausingFactor; ThetaArray; AngularDistribution; @@ -85,6 +84,9 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable DebugMode; DoSave; + SaveDirectory; + + Results; end properties (SetAccess = private, GetAccess = public) @@ -128,6 +130,8 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable @islogical); addParameter(p, 'SaveData', false,... @islogical); + addParameter(p, 'SaveDirectory', pwd,... + @ischar); p.parse(varargin{:}); @@ -144,7 +148,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable s.DebugMode = p.Results.DebugMode; s.DoSave = p.Results.SaveData; - + s.SaveDirectory = p.Results.SaveDirectory; s.reinitializeSimulator(); @@ -519,12 +523,6 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable function ret = get.VelocityCutoff(this) ret = this.VelocityCutoff; end - function set.CaptureFraction(this,val) - this.CaptureFraction = val; - end - function ret = get.CaptureFraction(this) - ret = this.CaptureFraction; - end function set.ClausingFactor(this,val) this.ClausingFactor = val; end @@ -570,6 +568,20 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable function ret = get.DoSave(this) ret = this.DoSave; 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 methods diff --git a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m index 75286c6..0472c67 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m +++ b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m @@ -11,7 +11,7 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector) Detuning = this.PushBeamDetuning - DopplerShift; 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 a_scatter = this.accelerationDueToSpontaneousEmissionProcess(SaturationIntensity, SaturationIntensity, Detuning, this.PushBeamLinewidth, this.PushBeamWaveVector); diff --git a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m index d1ffd2b..2f02c09 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m +++ b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m @@ -2,7 +2,7 @@ function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationInten Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1]; 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); if NumberOfScatteringEvents > 0 diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m b/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m index aed89df..98b4a67 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m @@ -20,6 +20,5 @@ function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector) break; end end - clear Index end diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m index 398e502..84cd469 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m @@ -1,4 +1,4 @@ -function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFraction, ClausingFactor, FinalDynamicalQuantities) +function [LoadingRate, StandardError] = calculateLoadingRate(this, FinalDynamicalQuantities) NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms); AutocorrelationFunction = zeros(1, this.NumberOfAtoms); @@ -23,41 +23,58 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti end for i = 1:this.NumberOfAtoms-1 - MeanLoadingRate = 0; + MeanTrappingEfficiency = 0; MeanLoadingRateShifted = 0; for j = 1:this.NumberOfAtoms-i - MeanLoadingRate = MeanLoadingRate + NumberOfLoadedAtoms(j)/j; + 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 * MeanLoadingRate * MeanLoadingRateShifted)); + AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanTrappingEfficiency * MeanLoadingRateShifted)); 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); - x = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms); - [FitObject,~] = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100); - n = FitObject.n; % n is the autocorrelation factor - MeanLoadingRate = 0; - NumberOfBins = floor(this.NumberOfAtoms/(2*n+1)); - LoadingRateError = zeros(1,NumberOfBins); - BinNumberLimit = min(NumberOfBins-1,5); - for i = 1:NumberOfBins-BinNumberLimit - LoadingRateError(i) = NumberOfLoadedAtoms(this.NumberOfAtoms-ceil((i-1)*(2*n+1))) / ... - (this.NumberOfAtoms-ceil((i-1)*(2*n+1))); - MeanLoadingRate = MeanLoadingRate + LoadingRateError(i); + 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 - MeanLoadingRate = MeanLoadingRate /(NumberOfBins-BinNumberLimit); - - StandardError = 0; - for i = 1:NumberOfBins-BinNumberLimit - StandardError = StandardError + (MeanLoadingRate-LoadingRateError(i))^2; + 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(StandardError/(NumberOfBins-BinNumberLimit)); - - LoadingRate = (MeanLoadingRate * this.calculateFreeMolecularRegimeFlux() * CaptureFraction * ClausingFactor); + + StandardError = sqrt((((NumberOfBlocks-BlockNumberLimit) - 1) / (NumberOfBlocks-BlockNumberLimit)) * Variance); else LoadingRate = 0; diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m index 7f0c6d0..d10ffa8 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m @@ -6,9 +6,7 @@ function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVecto %One side of parallelogram = PositionVectorFromWaveVectorOrigin %Base = Wavevector %Area = One side of parallelogram X Base - DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector); % Slow - %DistanceBetweenAtomAndLaserBeamAxis = norm((WaveVector*WaveVector')*PositionVectorFromWaveVectorOrigin-(WaveVector*PositionVectorFromWaveVectorOrigin')*WaveVector)./ ... - % (WaveVector(1)^2+WaveVector(2)^2+WaveVector(3)^2); % Faster + DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector); if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2); diff --git a/MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m b/MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m new file mode 100644 index 0000000..1deaf23 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m @@ -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 \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m b/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m index 243c0b2..a074cf0 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m +++ b/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m @@ -1,4 +1,4 @@ -function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ... +function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ... SecondParameterName, SecondParameterArray, varargin) p = inputParser; @@ -14,7 +14,10 @@ function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstPa addParameter(p, 'Order', 1, @(x) assert(isnumeric(x) && isscalar(x) && (x > 0) && (x < 3))); addParameter(p, 'RelatedParameterName', 'none', @ischar); addParameter(p, 'RelatedParameterArray', length(FirstParameterArray), @isvector); - + + addParameter(p, 'ChangeInitialConditions', false, @islogical); + addParameter(p, 'ParameterNameArray', {}, @iscell); + addParameter(p, 'ParameterValueArray', {}, @iscell); p.parse(this, FirstParameterName, FirstParameterArray, ... SecondParameterName, SecondParameterArray, varargin{:}) @@ -27,9 +30,14 @@ function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstPa Order = p.Results.Order; RelatedParameterName = p.Results.RelatedParameterName; RelatedParameterArray = p.Results.RelatedParameterArray; + ChangeInitialConditions = p.Results.ChangeInitialConditions; + ParameterNameArray = p.Results.ParameterNameArray; + ParameterValueArray = p.Results.ParameterValueArray; - NumberOfPointsForFirstParam = length(FirstParameterArray); - NumberOfPointsForSecondParam = length(SecondParameterArray); + NumberOfPointsForFirstParam = length(FirstParameterArray); + NumberOfPointsForSecondParam = length(SecondParameterArray); + LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); + StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); for i=1:NumberOfPointsForFirstParam eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i))); @@ -41,12 +49,30 @@ function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstPa eval(sprintf('OptionsStruct.%s = %d;', SecondParameterName, SecondParameterArray(j))); if ChangeRelatedParameter && Order == 2 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); this.setInitialConditions(options{:}); tic - [LoadingRate, ~] = this.runSimulation(); + [LoadingRate, StandardError] = this.runSimulation(); LoadingRateArray(i,j) = LoadingRate; + StandardErrorArray(i,j) = StandardError; 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 \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m b/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m index 809220a..49135f7 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m +++ b/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m @@ -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.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. - - 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 % Compute the angular distribution of the atomic beam @@ -58,5 +51,5 @@ function [LoadingRate, StandardError] = runSimulation(this) end clear Index %% Calculate the Loading Rate - [LoadingRate, StandardError] = this.calculateLoadingRate(this.CaptureFraction, this.ClausingFactor, FinalDynamicalQuantities); + [LoadingRate, StandardError] = this.calculateLoadingRate(FinalDynamicalQuantities); end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m b/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m index f5edc79..11cb451 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m +++ b/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m @@ -66,6 +66,7 @@ function setInitialConditions(this, varargin) this.CoolingBeamWaveVector = this.BlueWaveVector; this.CoolingBeamDetuning = this.BlueDetuning; this.CoolingBeamRadius = this.BlueBeamRadius; + this.CoolingBeamWaist = this.BlueBeamWaist; this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity; this.SidebandBeamRadius = this.BlueBeamRadius; this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity; diff --git a/MOT Capture Process Simulation/test_MOTSimulator.m b/MOT Capture Process Simulation/test_MOTSimulator.m index e0af22c..dcef493 100644 --- a/MOT Capture Process Simulation/test_MOTSimulator.m +++ b/MOT Capture Process Simulation/test_MOTSimulator.m @@ -6,13 +6,14 @@ OptionsStruct = struct; OptionsStruct.SimulationMode = '2D'; OptionsStruct.TimeStep = 50e-06; % in s OptionsStruct.SimulationTime = 4e-03; % in s -OptionsStruct.SpontaneousEmission = false; -OptionsStruct.Sideband = false; +OptionsStruct.SpontaneousEmission = true; +OptionsStruct.Sideband = true; OptionsStruct.ZeemanSlowerBeam = false; OptionsStruct.Gravity = true; OptionsStruct.AtomicBeamCollision = true; 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); Simulator = MOTSimulator(options{:}); @@ -52,7 +53,7 @@ Plotting.visualizeMagneticField(Simulator, XAxisRange, YAxisRange, ZAxisRange) %% - Plot MFP & VP for different temperatures TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius -Plotting.plotMeanFreePathAndVapourPressure(TemperatureinCelsius) +Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius) %% - Plot the Free Molecular Flux for different temperatures Temperature = [900, 950, 1000, 1050, 1100]; % Temperature' @@ -82,45 +83,54 @@ IncidentAtomPosition = 0; Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition) %% - 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 -NumberOfPointsForSecondParam = 2; -LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); +NumberOfPointsForParam = 20; %iterations of the simulation +% Scan Cooling Beam Power +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 -DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam); -PowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam); - -OptionsStruct.NumberOfAtoms = 5000; +OptionsStruct = struct; +OptionsStruct.ChangeInitialConditions = true; +OptionsStruct.ParameterNameArray = {'NumberOfAtoms'}; +OptionsStruct.ParameterValueArray = {10000}; +options = Helper.convertstruct2cell(OptionsStruct); tStart = tic; -for i=1:NumberOfPointsForFirstParam - 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 +[LoadingRateArray, ~] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:}); tEnd = toc(tStart); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); -%% TWO-PARAMETER SCAN FUNCTION +clear OptionsStruct -NumberOfPointsForFirstParam = 2; %iterations of the simulation -NumberOfPointsForSecondParam = 2; -Simulator.NumberOfAtoms = 5000; +% - Plot results + +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 -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; 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.RelatedParameterName = 'BluePower'; OptionsStruct.RelatedParameterArray = BluePowerArray; +OptionsStruct.ChangeInitialConditions = true; +OptionsStruct.ParameterNameArray = {'NumberOfAtoms'}; +OptionsStruct.ParameterValueArray = {10000}; options = Helper.convertstruct2cell(OptionsStruct); tStart = tic; -LoadingRateArray = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:}); +[LoadingRateArray, StandardErrorArray] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:}); tEnd = toc(tStart); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); clear OptionsStruct -%% - Plot results -FirstParameterArray = DetuningArray * Helper.PhysicsConstants.BlueLinewidth; -SecondParameterArray = (Simulator.TotalPower - (PowerArray * Simulator.TotalPower)); +% - Plot results + +FirstParameterArray = DetuningArray; +SecondParameterArray = SidebandPowerArray; QuantityOfInterestArray = LoadingRateArray; OptionsStruct = struct; @@ -148,7 +162,8 @@ OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.B OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)'; OptionsStruct.RescalingFactorForSecondParameter = 1000; 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); options = Helper.convertstruct2cell(OptionsStruct);