diff --git a/+Plotter/plotResultForThreeParameterScan.m b/+Plotter/plotResultForThreeParameterScan.m new file mode 100644 index 0000000..6867c48 --- /dev/null +++ b/+Plotter/plotResultForThreeParameterScan.m @@ -0,0 +1,82 @@ +function plotResultForThreeParameterScan(XParameter, YParameter, DeltaParameter, ZQuantity, varargin) + + p = inputParser; + p.KeepUnmatched = true; + + addRequired(p, 'FirstParameterArray', @isvector) + addRequired(p, 'SecondParameterArray', @isvector) + addRequired(p, 'ThirdParameterArray', @isvector) + addRequired(p, 'QuantityOfInterestArray', @ismatrix) + + addParameter(p, 'RescalingFactorForFirstParameter', 1, @isscalar) + addParameter(p, 'XLabelString', 'X parameter', @ischar) + addParameter(p, 'RescalingFactorForSecondParameter', 1, @isscalar) + addParameter(p, 'YLabelString', 'Y parameter', @ischar) + addParameter(p, 'RescalingFactorForThirdParameter', 1, @isscalar) + addParameter(p, 'RescalingFactorForQuantityOfInterest', 1, @isscalar) + addParameter(p, 'ZLabelString', 'Z parameter', @ischar) + addParameter(p, 'PlotTitleString', '', @ischar) + addParameter(p, 'FigureTitleString', 'Third-Parameter Scan', @ischar) + + p.parse(XParameter, YParameter, DeltaParameter, ZQuantity, varargin{:}) + + XParameter = p.Results.FirstParameterArray; + RescalingFactorForXParameter = p.Results.RescalingFactorForFirstParameter; + XLabelString = p.Results.XLabelString; + YParameter = p.Results.SecondParameterArray; + RescalingFactorForYParameter = p.Results.RescalingFactorForSecondParameter; + YLabelString = p.Results.YLabelString; + DeltaParameter = p.Results.ThirdParameterArray; + RescalingFactorForDeltaParameter = p.Results.RescalingFactorForThirdParameter; + ZQuantity = p.Results.QuantityOfInterestArray; + RescalingFactorForZQuantity = p.Results.RescalingFactorForQuantityOfInterest; + ZLabelString = p.Results.ZLabelString; + PlotTitleString = p.Results.PlotTitleString; + FigureTitleString = p.Results.FigureTitleString; + + f_h = Helper.getFigureByTag('Three-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 = 'Three-Parameter Scan'; + f_h.Units = 'pixels'; + + set(0,'units','pixels'); + screensize = get(0,'ScreenSize'); + f_h.Position = [[screensize(3)/60 screensize(4)/10] 1530 870]; + + RescaledXParameter = XParameter .* RescalingFactorForXParameter; + RescaledYParameter = YParameter .* RescalingFactorForYParameter; + RescaledDeltaParameter = DeltaParameter .* RescalingFactorForDeltaParameter; + + tiledlayout(ceil(length(DeltaParameter)/3), 3) + + for i = 1:length(DeltaParameter) + nexttile + RescaledZQuantity = ZQuantity{i} .* RescalingFactorForZQuantity; + imagesc(RescaledXParameter, RescaledYParameter, RescaledZQuantity(:,:)'); + set(gca,'YDir','normal'); + hXLabel = xlabel(XLabelString); + hYLabel = ylabel(YLabelString); + hPlotLabel = title(sprintf(PlotTitleString, RescaledDeltaParameter(i))); + set([hXLabel, hYLabel, hPlotLabel] , ... + 'FontSize' , 14); + end + + caxis([min(min(min(RescaledZQuantity))) max(max(max(RescaledZQuantity)))]); + + shading flat; + c = colorbar; + c.Label.String= ZLabelString; + c.Label.FontSize = 14; + c.Layout.Tile = 'east'; + + hTitle = sgtitle(FigureTitleString); + + set( hTitle , ... + 'FontSize' , 18 ); + + Helper.bringFiguresWithTagInForeground(); +end \ No newline at end of file diff --git a/+Simulator/+Scan/doOneParameter.m b/+Simulator/+Scan/doOneParameter.m new file mode 100644 index 0000000..5ac8b56 --- /dev/null +++ b/+Simulator/+Scan/doOneParameter.m @@ -0,0 +1,27 @@ +function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOneParameter(ovenObj, MOTobj, BeamName, BeamParameter, ParameterArray) + + Beams = MOTobj.Beams; + Beam = Beams{cellfun(@(x) strcmpi(x.Alias, BeamName), Beams)}; + NumberOfPointsForParam = length(ParameterArray); + LoadingRateArray = zeros(1,NumberOfPointsForParam); + StandardErrorArray = zeros(1,NumberOfPointsForParam); + ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2); + + for i=1:NumberOfPointsForParam + eval(sprintf('Beam.%s = %d;', BeamParameter, ParameterArray(i))); + [LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i,:)] = MOTobj.runSimulation(ovenObj); + end + + if MOTobj.DoSave + LoadingRate = struct; + LoadingRate.Values = LoadingRateArray; + LoadingRate.Errors = StandardErrorArray; + LoadingRate.CI = ConfidenceIntervalArray; + MOTobj.Results = LoadingRate; + SaveFolder = [MOTobj.SaveDirectory filesep 'Results']; + Filename = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')]; + eval([sprintf('%s_Object', Filename) ' = MOTobj;']); + mkdir(SaveFolder); + save([SaveFolder filesep Filename], sprintf('%s_Object', Filename)); + end +end \ No newline at end of file diff --git a/+Simulator/+Scan/doThreeParameters.m b/+Simulator/+Scan/doThreeParameters.m new file mode 100644 index 0000000..4afdcf3 --- /dev/null +++ b/+Simulator/+Scan/doThreeParameters.m @@ -0,0 +1,22 @@ +function LoadingRateArray = doThreeParameters(ovenObj, MOTobj, BeamName, FirstBeamParameter, FirstParameterArray, ... + SecondBeamParameter, SecondParameterArray, ThirdBeamParameter, ThirdParameterArray) + + NumberOfPointsForThirdParam = length(ThirdParameterArray); + LoadingRateArray = {}; + + for i=1:NumberOfPointsForThirdParam + eval(sprintf('MOTobj.%s = %d;', ThirdBeamParameter, ThirdParameterArray(i))); + LoadingRateArray{end+1} = Simulator.Scan.doTwoParameters(ovenObj, MOTobj, BeamName, FirstBeamParameter, FirstParameterArray, SecondBeamParameter, SecondParameterArray); + end + + if MOTobj.DoSave + LoadingRate = struct; + LoadingRate.Values = LoadingRateArray; + MOTobj.Results = LoadingRate; + SaveFolder = [MOTobj.SaveDirectory filesep 'Results']; + Filename = ['ThreeParameterScan_' datestr(now,'yyyymmdd_HHMM')]; + eval([sprintf('%s_Object', Filename) ' = MOTobj;']); + mkdir(SaveFolder); + save([SaveFolder filesep Filename], sprintf('%s_Object', Filename)); + end +end \ No newline at end of file diff --git a/+Simulator/+Scan/doTwoParameters.m b/+Simulator/+Scan/doTwoParameters.m new file mode 100644 index 0000000..092ca13 --- /dev/null +++ b/+Simulator/+Scan/doTwoParameters.m @@ -0,0 +1,31 @@ +function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwoParameters(ovenObj, MOTobj, BeamName, FirstBeamParameter, ... + FirstParameterArray, SecondBeamParameter, SecondParameterArray) + Beams = MOTobj.Beams; + Beam = Beams{cellfun(@(x) strcmpi(x.Alias, BeamName), Beams)}; + NumberOfPointsForFirstParam = length(FirstParameterArray); + NumberOfPointsForSecondParam = length(SecondParameterArray); + LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); + StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); + ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2); + + for i=1:NumberOfPointsForFirstParam + eval(sprintf('Beam.%s = %d;', FirstBeamParameter, FirstParameterArray(i))); + for j=1:NumberOfPointsForSecondParam + eval(sprintf('Beam.%s = %d;', SecondBeamParameter, SecondParameterArray(j))); + [LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOTobj.runSimulation(ovenObj); + end + end + + if MOTobj.DoSave + LoadingRate = struct; + LoadingRate.Values = LoadingRateArray; + LoadingRate.Errors = StandardErrorArray; + LoadingRate.CI = ConfidenceIntervalArray; + MOTobj.Results = LoadingRate; + SaveFolder = [MOTobj.SaveDirectory filesep 'Results']; + Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')]; + eval([sprintf('%s_Object', Filename) ' = MOTobj;']); + mkdir(SaveFolder); + save([SaveFolder filesep Filename], sprintf('%s_Object', Filename)); + end +end \ No newline at end of file diff --git a/test_MOTCaptureProcessSimulation.m b/test_MOTCaptureProcessSimulation.m index a74b82a..624f967 100644 --- a/test_MOTCaptureProcessSimulation.m +++ b/test_MOTCaptureProcessSimulation.m @@ -7,7 +7,6 @@ % - Automatically creates Beams objects OptionsStruct = struct; OptionsStruct.ErrorEstimationMethod = 'bootstrap'; % 'jackknife' | 'bootstrap' -OptionsStruct.NumberOfAtoms = 10000; OptionsStruct.TimeStep = 50e-06; % in s OptionsStruct.SimulationTime = 4e-03; % in s OptionsStruct.SpontaneousEmission = true; @@ -15,8 +14,8 @@ OptionsStruct.Sideband = false; OptionsStruct.PushBeam = true; OptionsStruct.Gravity = true; OptionsStruct.BackgroundCollision = true; -OptionsStruct.SaveData = false; -OptionsStruct.SaveDirectory = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation'; +OptionsStruct.SaveData = true; +% OptionsStruct.SaveDirectory = ''; options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct @@ -70,7 +69,13 @@ Plotter.plotAngularDistributionForDifferentBeta(Oven, Beta) Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot! %% - Plot Phase Space with Acceleration Field -MOT2D.SidebandBeam = false; +MOT2D.SidebandBeam = true; +CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; +CoolingBeam.Power = 0.2; +CoolingBeam.Detuning = -1.67*Helper.PhysicsConstants.BlueLinewidth; +SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)}; +SidebandBeam.Power = 0.2; +SidebandBeam.Detuning = -3.35*Helper.PhysicsConstants.BlueLinewidth; MOT2D.NumberOfAtoms = 50; MinimumVelocity = 0; MaximumVelocity = 150; @@ -93,35 +98,18 @@ Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirect %% - Scan parameters: One-Parameter Scan -MOT2D.NumberOfAtoms = 5000; -MOT2D.TotalPower = 0.4; -CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; - +MOT2D.NumberOfAtoms = 5000; +MOT2D.TotalPower = 0.4; NumberOfPointsForFirstParam = 5; %iterations of the simulation -% Scan Cooling Beam Power -PowerArray = linspace(0.1, 1.0, NumberOfPointsForFirstParam) * MOT2D.TotalPower; -% Scan Cooling Beam Detuning -% DetuningArray = linspace(-0.5,-10, NumberOfPointsForParam) * Helper.PhysicsConstants.BlueLinewidth; - -LoadingRateArray = zeros(1,NumberOfPointsForFirstParam); -StandardErrorArray = zeros(1,NumberOfPointsForFirstParam); -ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, 2); +ParameterArray = linspace(0.1, 1.0, NumberOfPointsForFirstParam) * MOT2D.TotalPower; tStart = tic; -for i=1:NumberOfPointsForFirstParam - CoolingBeam.Power = PowerArray(i); - [LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i, :)] = MOT2D.runSimulation(Oven); -end +[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.Scan.doOneParameter(Oven, MOT2D, 'Blue', 'Power', ParameterArray); tEnd = toc(tStart); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); -clear OptionsStruct - % - Plot results -ParameterArray = PowerArray; -QuantityOfInterestArray = LoadingRateArray; - OptionsStruct = struct; OptionsStruct.RescalingFactorForParameter = 1000; OptionsStruct.XLabelString = 'Cooling Beam Power (mW)'; @@ -136,87 +124,112 @@ OptionsStruct.TitleString = sprintf('Magnetic Gradient options = Helper.convertstruct2cell(OptionsStruct); -Plotter.plotResultForOneParameterScan(ParameterArray, QuantityOfInterestArray, options{:}) +Plotter.plotResultForOneParameterScan(ParameterArray, LoadingRateArray, options{:}) clear OptionsStruct -if MOT2D.DoSave - LoadingRate = struct; - LoadingRate.Values = LoadingRateArray; - LoadingRate.Errors = StandardErrorArray; - LoadingRate.CI = ConfidenceIntervalArray; - MOT2D.Results = LoadingRate; - SaveFolder = [MOT2D.SaveDirectory filesep 'Results']; - Filename = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')]; - eval([sprintf('%s_Object', Filename) ' = MOT2D;']); - mkdir(SaveFolder); - save([SaveFolder filesep Filename], sprintf('%s_Object', Filename)); -end - %% - Scan parameters: Two-Parameter Scan -MOT2D.NumberOfAtoms = 50; -MOT2D.TotalPower = 0.6; -MOT2D.SidebandBeam = false; -SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)}; +% COOLING BEAM POWER VS DETUNING +MOT2D.NumberOfAtoms = 5000; +MOT2D.TotalPower = 0.6; NumberOfPointsForFirstParam = 10; %iterations of the simulation NumberOfPointsForSecondParam = 10; - -% Scan Sideband Detuning and Power Ratio -DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; -% SidebandPowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam) * MOT2D.TotalPower; -% BluePowerArray = MOT2D.TotalPower - SidebandPowerArray; -BluePowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam) * MOT2D.TotalPower; - -LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); -StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); -ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2); +FirstParameterArray = linspace(-0.5, -2.5, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; +SecondParameterArray = linspace(0.3, 1.0, NumberOfPointsForSecondParam) * MOT2D.TotalPower; tStart = tic; -for i = 1:NumberOfPointsForFirstParam - SidebandBeam.Detuning = DetuningArray(i); - for j = 1:NumberOfPointsForSecondParam - SidebandBeam.Power = SidebandPowerArray(j); - CoolingBeam.Power = BluePowerArray(j); - [LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOT2D.runSimulation(Oven); - end -end +[LoadingRateArray, ~, ~] = Simulator.Scan.doTwoParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, 'Power', SecondParameterArray); tEnd = toc(tStart); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); -clear OptionsStruct - % - Plot results -FirstParameterArray = DetuningArray; -SecondParameterArray = SidebandPowerArray; -QuantityOfInterestArray = LoadingRateArray; - OptionsStruct = struct; OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1; -OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)'; +OptionsStruct.XLabelString = 'Cooling Beam Detuning (\Delta/\Gamma)'; OptionsStruct.RescalingFactorForSecondParameter = 1000; -OptionsStruct.YLabelString = 'Sideband Power (mW)'; -OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-11; -OptionsStruct.ZLabelString = 'Loading rate (x 10^{11} atoms/s)'; +OptionsStruct.YLabelString = 'Cooling Beam Power (mW)'; +OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-9; +OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)'; OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100); options = Helper.convertstruct2cell(OptionsStruct); -Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, QuantityOfInterestArray, options{:}) +Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, LoadingRateArray, options{:}) clear OptionsStruct -if MOT2D.DoSave - LoadingRate = struct; - LoadingRate.Values = LoadingRateArray; - LoadingRate.Errors = StandardErrorArray; - LoadingRate.CI = ConfidenceIntervalArray; - MOT2D.Results = LoadingRate; - SaveFolder = [MOT2D.SaveDirectory filesep 'Results']; - Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')]; - eval([sprintf('%s_Object', Filename) ' = MOT2D;']); - mkdir(SaveFolder); - save([SaveFolder filesep Filename], sprintf('%s_Object', Filename)); -end \ No newline at end of file +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% COOLING BEAM WAIST VS DETUNING + +MOT2D.NumberOfAtoms = 5000; +CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; +CoolingBeam.Power = 0.4; +NumberOfPointsForFirstParam = 10; %iterations of the simulation +NumberOfPointsForSecondParam = 10; +FirstParameterArray = linspace(-0.5, -2.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; +SecondParameterArray = linspace(10, 25, NumberOfPointsForSecondParam) * 1e-03; + +tStart = tic; +[LoadingRateArray, ~, ~] = Simulator.Scan.doTwoParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, 'Waist', SecondParameterArray); +tEnd = toc(tStart); +fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); + +% - Plot results + +OptionsStruct = struct; +OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1; +OptionsStruct.XLabelString = 'Cooling Beam Detuning (\Delta/\Gamma)'; +OptionsStruct.RescalingFactorForSecondParameter = 1000; +OptionsStruct.YLabelString = 'Cooling Beam Waist (mm)'; +OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-9; +OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)'; +OptionsStruct.TitleString = sprintf('Cooling Beam Power = %d (mW); Magnetic Gradient = %.0f (G/cm)', CoolingBeam.Power*1000, MOT2D.MagneticGradient * 100); + +options = Helper.convertstruct2cell(OptionsStruct); + +Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, LoadingRateArray, options{:}) + +clear OptionsStruct + +%% - Scan parameters: Three-Parameter Scan + +% COOLING BEAM WAIST VS DETUNING FOR DIFFERENT MAGNETIC FIELD GRADIENTS + +MOT2D.NumberOfAtoms = 5000; +CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; +CoolingBeam.Power = 0.4; +NumberOfPointsForFirstParam = 10; %iterations of the simulation +NumberOfPointsForSecondParam = 10; +NumberOfPointsForThirdParam = 6; +FirstParameterArray = linspace(-0.5, -2.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; +SecondParameterArray = linspace(10, 25, NumberOfPointsForSecondParam) * 1e-03; +ThirdParameterArray = linspace(30, 50, NumberOfPointsForThirdParam) * 1e-02; + +tStart = tic; +LoadingRateArray = Simulator.Scan.doThreeParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, ... + 'Waist', SecondParameterArray, ... + 'MagneticGradient', ThirdParameterArray); +tEnd = toc(tStart); +fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); + +% - Plot results + +OptionsStruct = struct; +OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1; +OptionsStruct.XLabelString = 'Cooling Beam Detuning (\Delta/\Gamma)'; +OptionsStruct.RescalingFactorForSecondParameter = 1000; +OptionsStruct.YLabelString = 'Cooling Beam Waist (mm)'; +OptionsStruct.RescalingFactorForThirdParameter = 100; +OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-9; +OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)'; +OptionsStruct.PlotTitleString = 'Magnetic Gradient = %.0f (G/cm)'; +OptionsStruct.FigureTitleString = sprintf('Oven-2DMOT Distance = %.1f (mm); Cooling Beam Power = %d (mW)', Oven.OvenDistance * 1000, CoolingBeam.Power*1000); + +options = Helper.convertstruct2cell(OptionsStruct); + +Plotter.plotResultForThreeParameterScan(FirstParameterArray, SecondParameterArray, ThirdParameterArray, LoadingRateArray, options{:}) + +clear OptionsStruct \ No newline at end of file