From add5d02ae3cbe349f2cacf72065bd8fc7f78e9b6 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sat, 3 Jul 2021 10:19:27 +0200 Subject: [PATCH] Added new plotting routines that allow for checks, minor cosmetic changes to other scripts. --- .../+Helper/PhysicsConstants.m | 1 + .../+Helper/bringFiguresWithTagInForeground.m | 15 ++++ .../+Helper/getFigureByTag.m | 2 +- .../plotAngularDistributionForDifferentBeta.m | 76 +++++++++++++++++ .../+Plotting/plotCaptureVelocityVsAngle.m | 37 ++++++++ .../+Plotting/plotFreeMolecularFluxVsTemp.m | 54 ++++++++++++ .../plotMeanFreePathAndVapourPressureVsTemp.m | 43 ++++++++++ .../plotPhaseSpaceWithAccelerationField.m | 81 ++++++++++++++++++ .../plotPositionAndVelocitySampling.m | 8 +- .../+Plotting/plotResultForTwoParameterScan.m | 2 + .../+Plotting/visualizeMagneticField.m | 2 + .../@MOTSimulator/MOTSimulator.m | 41 ++++++--- .../@MOTSimulator/calculateCaptureVelocity.m | 1 + .../calculateFreeMolecularRegimeFlux.m | 7 +- .../@MOTSimulator/calculateLoadingRate.m | 18 ++-- .../calculateLocalSaturationIntensity.m | 6 +- .../calculateTotalAcceleration.m | 2 +- .../@MOTSimulator/doTwoParameterScan.m | 52 ++++++++++++ .../@MOTSimulator/initialVelocitySampling.m | 24 +++--- .../@MOTSimulator/reinitializeSimulator.m | 9 +- .../@MOTSimulator/runSimulation.m | 12 +-- .../@MOTSimulator/setInitialConditions.m | 85 ++++++++++++++----- .../@MOTSimulator/solver.m | 4 +- .../test_MOTSimulator.m | 78 ++++++++++++++--- 24 files changed, 572 insertions(+), 88 deletions(-) create mode 100644 MOT Capture Process Simulation/+Helper/bringFiguresWithTagInForeground.m create mode 100644 MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m create mode 100644 MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m create mode 100644 MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m create mode 100644 MOT Capture Process Simulation/+Plotting/plotMeanFreePathAndVapourPressureVsTemp.m create mode 100644 MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m diff --git a/MOT Capture Process Simulation/+Helper/PhysicsConstants.m b/MOT Capture Process Simulation/+Helper/PhysicsConstants.m index 65efdae..8262415 100644 --- a/MOT Capture Process Simulation/+Helper/PhysicsConstants.m +++ b/MOT Capture Process Simulation/+Helper/PhysicsConstants.m @@ -24,6 +24,7 @@ classdef PhysicsConstants < handle % Dy specific constants Dy164Mass = 163.929174751*1.66053878283E-27; + Dy164IsotopicAbundance = 0.2826; BlueWavelength = 421.291e-9; BlueLandegFactor = 1.22; BlueLifetime = 4.94e-9; diff --git a/MOT Capture Process Simulation/+Helper/bringFiguresWithTagInForeground.m b/MOT Capture Process Simulation/+Helper/bringFiguresWithTagInForeground.m new file mode 100644 index 0000000..c58a117 --- /dev/null +++ b/MOT Capture Process Simulation/+Helper/bringFiguresWithTagInForeground.m @@ -0,0 +1,15 @@ +function output = bringFiguresWithTagInForeground() + +figure_handles = findobj('type','figure'); + +for idx = 1:length(figure_handles) + if ~isempty(figure_handles(idx).Tag) + figure(figure_handles(idx)); + end +end + +if nargout > 0 + output = figure_handles; +end + +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Helper/getFigureByTag.m b/MOT Capture Process Simulation/+Helper/getFigureByTag.m index 7932a97..8fc6bf9 100644 --- a/MOT Capture Process Simulation/+Helper/getFigureByTag.m +++ b/MOT Capture Process Simulation/+Helper/getFigureByTag.m @@ -183,7 +183,7 @@ function copyToClipboard(~,~) end else pos = fig_h.Position; - Helper.screencapture(fig_h,[],'clipboard','position',[1,1,pos(3)-2,pos(4)]); + Helper.screencapture(fig_h,[],'clipboard','position',[7,7,pos(3)-2,pos(4)]); end end diff --git a/MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m b/MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m new file mode 100644 index 0000000..5b9f7db --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m @@ -0,0 +1,76 @@ +function plotAngularDistributionForDifferentBeta(obj, Beta) + + f_h = Helper.getFigureByTag('AngDistForBeta'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Beta dependence'; + 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]; + + hold on + + obj.reinitializeSimulator(); + + SimulationBeta = obj.Beta; + + if ~ismember(SimulationBeta, Beta) + theta = linspace(0.0001,pi/2,1000); + J = zeros(1,length(theta)); + for j=1:length(theta) + J(j) = obj.angularDistributionFunction(theta(j)); + end + Norm = 0; + for j=1:length(J) + Norm = Norm+J(j)*sin(theta(j))*(theta(2)-theta(1)); + end + J = J ./Norm*2; + J = J ./max(J); + plot([-flip(theta),theta], [flip(J),J],'DisplayName', ['\beta = ',num2str(SimulationBeta, '%.3f')], 'Linewidth',1.5) + end + + for i=1:length(Beta) + + theta = linspace(0.0001,pi/2,1000); + J = zeros(1,length(theta)); + obj.Beta = Beta(i); + obj.NozzleLength = (2 * obj.NozzleRadius) / Beta(i); + + for j=1:length(theta) + J(j) = obj.angularDistributionFunction(theta(j)); + end + + Norm = 0; + for j=1:length(J) + Norm = Norm+J(j)*sin(theta(j))*(theta(2)-theta(1)); + end + J = J ./Norm*2; + J = J ./max(J); + + if Beta(i) ~= SimulationBeta + plot([-flip(theta),theta], [flip(J),J],'DisplayName',['\beta = ',num2str(Beta(i))], 'LineStyle', '--', 'Linewidth',1.5) + else + plot([-flip(theta),theta], [flip(J),J],'DisplayName',['\beta = ',num2str(Beta(i))], 'Linewidth',1.5) + end + end + + hold off + + leg = legend; + hXLabel = xlabel('\theta (rad)'); + hYLabel = ylabel('J(\theta)'); + hTitle = sgtitle('Angular Distribution (Free Molecular Flow)'); + + set([hXLabel, hYLabel, leg] , ... + 'FontSize' , 14 ); + set( hTitle , ... + 'FontSize' , 18 ); + + grid on + Helper.bringFiguresWithTagInForeground(); +end diff --git a/MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m b/MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m new file mode 100644 index 0000000..bfff80e --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m @@ -0,0 +1,37 @@ +function plotCaptureVelocityVsAngle(obj) + + theta = linspace(0, obj.MOTExitDivergence, 40); + CaptureVelocity = zeros(length(theta),3); + + for i=1:length(theta) + CaptureVelocity(i,:) = obj.calculateCaptureVelocity([-obj.OvenDistance,0,0], [cos(theta(i)),0,sin(theta(i))]); + end + + f_h = Helper.getFigureByTag('Capture Velocity'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Capture Velocity'; + 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]; + + plot(theta, sqrt(CaptureVelocity(:,1).^2+CaptureVelocity(:,2).^2+CaptureVelocity(:,3).^2), 'Linewidth', 1.5) + + hXLabel = xlabel('\theta (rad)'); + hYLabel = ylabel('Velocity (m/s)'); + hTitle = sgtitle('Capture Velocity of an atomic beam from (0,0,0)'); + + set([hXLabel, hYLabel] , ... + 'FontSize' , 14 ); + set( hTitle , ... + 'FontSize' , 18 ); + + grid on + Helper.bringFiguresWithTagInForeground(); + +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m b/MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m new file mode 100644 index 0000000..b2edac2 --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m @@ -0,0 +1,54 @@ +function plotFreeMolecularFluxVsTemp(obj, Temperature) + + f_h = Helper.getFigureByTag('Tube Flux'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Tube Flux'; + f_h.Units = 'pixels'; + + set(0,'units','pixels'); + screensize = get(0,'ScreenSize'); + f_h.Position = [[screensize(3)/4.5 screensize(4)/4.5] 920 750]; + + hold on + + for i=1:length(Temperature) + beta = linspace(0.01,0.5,200); + L = 2*obj.NozzleRadius./beta; + obj.OvenTemperature = Temperature(i); + flux = zeros(1,length(L)); + for j=1:length(L) + obj.NozzleLength = L(j); + flux(j) = obj.calculateFreeMolecularRegimeFlux(); + end + plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5) + end + set(gca,'yscale','log') + + obj.reinitializeSimulator(); + + xline(obj.Beta, 'k--','Linewidth', 0.5); + fmf = obj.calculateFreeMolecularRegimeFlux(); + yline(fmf, 'k--', 'Linewidth', 1.5); + textstring = [sprintf('%1.e',fmf) ' atoms/s for ' '$$ \beta $$ = ' num2str(obj.Beta, '%.2f') sprintf(' @ %.2f K', obj.OvenTemperatureinKelvin)]; + txt = text((obj.Beta + 0.05*obj.Beta), (max(fmf) + 0.2*fmf), textstring, 'Interpreter','latex'); + + hold off + + leg = legend('Location', 'southeast'); + leg.String = leg.String(1:end-2); % Remove xline and yline legend entries + hXLabel = xlabel('\beta'); + hYLabel = ylabel('Flux (atoms/s)'); + hTitle = sgtitle('Total Flux from a Tube (Free Molecular Flow)'); + + set([hXLabel, hYLabel, txt, leg] , ... + 'FontSize' , 14 ); + set( hTitle , ... + 'FontSize' , 18 ); + + grid on + Helper.bringFiguresWithTagInForeground(); +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Plotting/plotMeanFreePathAndVapourPressureVsTemp.m b/MOT Capture Process Simulation/+Plotting/plotMeanFreePathAndVapourPressureVsTemp.m new file mode 100644 index 0000000..e7b9551 --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotMeanFreePathAndVapourPressureVsTemp.m @@ -0,0 +1,43 @@ +function plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius) + + TemperatureinKelvin = TemperatureinCelsius + 273.15; + + Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+TemperatureinKelvin)); % Vapor Pressure Dysprosium for the given oven temperature + MeanFreepath = (Helper.PhysicsConstants.BoltzmannConstant*TemperatureinKelvin)./(sqrt(2) * ( pi * (2*281e-12)^2) * Dy164VapourPressure); %free path at above tempeture + + f_h = Helper.getFigureByTag('Dysprosium Gas Properties'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Dysprosium Gas Properties'; + 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]; + + yyaxis left + semilogy(TemperatureinCelsius, Dy164VapourPressure, 'Color', '#0071BB', 'Linewidth',1.5); + hYLabel_1 = ylabel('Vapor Pressure (mbar)'); + yyaxis right + semilogy(TemperatureinCelsius, MeanFreepath.*1000, 'Color', '#36B449', 'Linewidth',1.5) + hYLabel_2 = ylabel('Mean Free Path (mm)'); + + hXLabel = xlabel('Temperature (℃)'); + + ax = gca; + ax.YAxis(1).Color = '#0071BB'; + ax.YAxis(2).Color = '#36B449'; + + hTitle = sgtitle('^{164}Dy Gas'); + + set([hXLabel, hYLabel_1, hYLabel_2] , ... + 'FontSize' , 14 ); + set( hTitle , ... + 'FontSize' , 18 ); + + grid on + Helper.bringFiguresWithTagInForeground(); +end diff --git a/MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m b/MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m new file mode 100644 index 0000000..7edb2c7 --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m @@ -0,0 +1,81 @@ +function plotPhaseSpaceWithAccelerationField(obj, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition) + + f_h = Helper.getFigureByTag('Phase Space Plot'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Phase Space Plot'; + 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]; + + N = obj.NumberOfAtoms; + L = obj.OvenDistance * 2; + Theta = IncidentAtomDirection; + z = IncidentAtomPosition; + T = obj.SimulationTime; + tau = obj.TimeStep; + + [X,Y] = meshgrid(-L/2:L/NumberOfBins:L/2,-MaximumVelocity:2*MaximumVelocity/NumberOfBins:MaximumVelocity); + + a=zeros(NumberOfBins+1,NumberOfBins+1,3); + + obj.setInitialConditions(); + + for i=1:length(X) + for j=1:length(Y) + a(i,j,:) = obj.calculateTotalAcceleration([X(1,i), 0, z], [Y(j,1)*cos(Theta),0,Y(j,1)*sin(Theta)]); + end + end + for i=1:length(X) + for j=1:length(Y) + if isnan(a(i,j,1)) || isnan(a(i,j,2)) || isnan(a(i,j,3)) + a(i,j,1)=0; + a(i,j,2)=0; + a(i,j,3)=0; + end + end + end + + pcolor(X',Y',a(:,:,1)) + hold on + col=colorbar; + col.Label.String='Aceleration (m/s^2)'; + shading flat + + %------------------------------------------------------------------------- + + Y = linspace(0,MaximumVelocity,N); + Trajectories = zeros(length(Y),int64(T/tau),9); + for i=1:length(Y) + x =-L/2; + vx = Y(i)*cos(Theta); + vz = Y(i)*sin(Theta); + r = [x,0,z]; + v = [vx,0,vz]; + [Trajectories(i,:,:),~] = obj.solver(r, v); + end + + hold on + + for i=1:length(Y) + plot(Trajectories(i,:,1),Trajectories(i,:,4),'w','linewidth',1.3) + end + + hold off + + hXLabel = xlabel('Position: Along the x-axis (m)'); + hYLabel = ylabel('Velocity (m/s)'); + hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient)); + set([hXLabel, hYLabel] , ... + 'FontSize' , 14 ); + set( hTitle , ... + 'FontSize' , 18 ); + + Helper.bringFiguresWithTagInForeground(); +end + diff --git a/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m b/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m index 0d2692a..d67a91c 100644 --- a/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m +++ b/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m @@ -1,4 +1,4 @@ -function plotPositionAndVelocitySampling(Simulator, NumberOfBins) +function plotPositionAndVelocitySampling(obj, NumberOfBins) f_h = Helper.getFigureByTag('RejectionSampling'); set(groot,'CurrentFigure',f_h); @@ -13,8 +13,8 @@ function plotPositionAndVelocitySampling(Simulator, NumberOfBins) screensize = get(0,'ScreenSize'); f_h.Position = [[screensize(3)/7 screensize(4)/7] 1.357e+03 770]; - initialPositions = Simulator.InitialPositions; - initialVelocities = Simulator.InitialVelocities; + initialPositions = obj.InitialPositions; + initialVelocities = obj.InitialVelocities; subplot(3,2,1) histogram(initialPositions(:, 1)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component') @@ -53,4 +53,6 @@ function plotPositionAndVelocitySampling(Simulator, NumberOfBins) legend('FontSize', 14) sgtitle('Rejection sampling for initial distributions','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 38ca88e..91bafd6 100644 --- a/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m +++ b/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m @@ -62,4 +62,6 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg 'FontSize' , 14 ); set( hTitle , ... 'FontSize' , 18 ); + + Helper.bringFiguresWithTagInForeground(); end \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m b/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m index 78f895a..2af9d23 100644 --- a/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m +++ b/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m @@ -67,4 +67,6 @@ function visualizeMagneticField(obj, x_range, y_range, z_range) 'FontSize' , 14 ); set( hTitle , ... 'FontSize' , 18 ); + + Helper.bringFiguresWithTagInForeground(); end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m index 35b85ec..9cb189b 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m +++ b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m @@ -72,6 +72,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable VelocityCutoff; CaptureFraction; ClausingFactor; + ThetaArray; AngularDistribution; NormalizationConstantForAngularDistribution; @@ -80,13 +81,14 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable Sideband; ZeemanSlowerBeam; Gravity; + AtomicBeamCollision; DebugMode; DoSave; end properties (SetAccess = private, GetAccess = public) - InitialParameters + SimulationParameters end properties (Dependent, SetAccess = private) @@ -120,6 +122,8 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable @islogical); addParameter(p, 'Gravity', false,... @islogical); + addParameter(p, 'AtomicBeamCollision', true,... + @islogical); addParameter(p, 'DebugMode', false,... @islogical); addParameter(p, 'SaveData', false,... @@ -127,18 +131,21 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable p.parse(varargin{:}); - s.SimulationMode = p.Results.SimulationMode; + s.SimulationMode = p.Results.SimulationMode; - s.TimeStep = p.Results.TimeStep; - s.SimulationTime = p.Results.SimulationTime; + s.TimeStep = p.Results.TimeStep; + s.SimulationTime = p.Results.SimulationTime; s.SpontaneousEmission = p.Results.SpontaneousEmission; s.Sideband = p.Results.Sideband; s.ZeemanSlowerBeam = p.Results.ZeemanSlowerBeam; s.Gravity = p.Results.Gravity; - - s.DebugMode = p.Results.DebugMode; - s.DoSave = p.Results.SaveData; + s.AtomicBeamCollision = p.Results.AtomicBeamCollision; + + s.DebugMode = p.Results.DebugMode; + s.DoSave = p.Results.SaveData; + + s.reinitializeSimulator(); @@ -159,7 +166,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable ret = this.TimeStep; end function set.SimulationTime(this, val) - assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!'); +% assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!'); this.SimulationTime = val; end function ret = get.SimulationTime(this) @@ -530,13 +537,27 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable function ret = get.AngularDistribution(this) ret = this.AngularDistribution; end + function set.ThetaArray(this,val) + this.ThetaArray = val; + end + function ret = get.ThetaArray(this) + ret = this.ThetaArray; + end function set.NormalizationConstantForAngularDistribution(this,val) this.NormalizationConstantForAngularDistribution = val; end function ret = get.NormalizationConstantForAngularDistribution(this) ret = this.NormalizationConstantForAngularDistribution; end - + + function set.AtomicBeamCollision(this,val) + this.AtomicBeamCollision=val; + end + + function ret=get.AtomicBeamCollision(this) + ret=this.AtomicBeamCollision; + end + function set.DebugMode(this, val) this.DebugMode = val; end @@ -544,7 +565,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable ret = this.DebugMode; end function set.DoSave(this, val) - this.DebugMode = val; + this.DoSave = val; end function ret = get.DoSave(this) ret = this.DoSave; diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m b/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m index 837388e..aed89df 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m @@ -3,6 +3,7 @@ function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector) VelocityVector = VelocityVector./norm(VelocityVector); UpperLimit = 500; LowerLimit = 0; + this.AtomicBeamCollision = false; for Index = 1:500 InitialVelocity = 0.5 * (UpperLimit + LowerLimit) * VelocityVector; diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m b/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m index 7a32584..5307c3a 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m @@ -1,7 +1,8 @@ function ret = calculateFreeMolecularRegimeFlux(this) %This function calculate the total flux of atoms coming out from a tube %See Expected atomic flux section in Barbiero - Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature - Dy164DensityinOven = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin); - ret =1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2; + Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature + Dy164DensityinOven = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin); + ClausingFactorApproximation = (8 * this.NozzleRadius) / (3 * this.NozzleLength); + ret = Helper.PhysicsConstants.Dy164IsotopicAbundance * 1/4 * Dy164DensityinOven * this.AverageVelocity * ClausingFactorApproximation * pi * this.NozzleRadius.^2; end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m index 3c2d225..398e502 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m @@ -3,7 +3,7 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms); AutocorrelationFunction = zeros(1, this.NumberOfAtoms); - for i = 1:NumberOfLoadedAtoms + for i = 1:this.NumberOfAtoms FinalPosition = FinalDynamicalQuantities(i,1:3); DivergenceAngle = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2))); if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0) @@ -22,30 +22,30 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti end end - for i = 1:NumberOfLoadedAtoms-1 + for i = 1:this.NumberOfAtoms-1 MeanLoadingRate = 0; MeanLoadingRateShifted = 0; - for j = 1:NumberOfLoadedAtoms-i + for j = 1:this.NumberOfAtoms-i MeanLoadingRate = MeanLoadingRate + 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) = ((NumberOfLoadedAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((NumberOfLoadedAtoms-i)^-1 * MeanLoadingRate * MeanLoadingRateShifted)); + AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanLoadingRate * MeanLoadingRateShifted)); end if AutocorrelationFunction(1)~=0 % In case no atom loading AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1); - x = linspace(1, NumberOfLoadedAtoms, NumberOfLoadedAtoms); + 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(NumberOfLoadedAtoms/(2*n+1)); + NumberOfBins = floor(this.NumberOfAtoms/(2*n+1)); LoadingRateError = zeros(1,NumberOfBins); BinNumberLimit = min(NumberOfBins-1,5); for i = 1:NumberOfBins-BinNumberLimit - LoadingRateError(i) = NumberOfLoadedAtoms(NumberOfLoadedAtoms-ceil((i-1)*(2*n+1))) / ... - (NumberOfLoadedAtoms-ceil((i-1)*(2*n+1))); + LoadingRateError(i) = NumberOfLoadedAtoms(this.NumberOfAtoms-ceil((i-1)*(2*n+1))) / ... + (this.NumberOfAtoms-ceil((i-1)*(2*n+1))); MeanLoadingRate = MeanLoadingRate + LoadingRateError(i); end @@ -57,7 +57,7 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti end StandardError = sqrt(StandardError/(NumberOfBins-BinNumberLimit)); - LoadingRate = (MeanLoadingRate * this.FreeMolecularRegimeFlux() * CaptureFraction * ClausingFactor); + LoadingRate = (MeanLoadingRate * this.calculateFreeMolecularRegimeFlux() * CaptureFraction * ClausingFactor); else LoadingRate = 0; diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m index 5472d29..7f0c6d0 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m @@ -6,9 +6,9 @@ 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); % Slow + %DistanceBetweenAtomAndLaserBeamAxis = norm((WaveVector*WaveVector')*PositionVectorFromWaveVectorOrigin-(WaveVector*PositionVectorFromWaveVectorOrigin')*WaveVector)./ ... + % (WaveVector(1)^2+WaveVector(2)^2+WaveVector(3)^2); % Faster if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2); diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m b/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m index b7da57c..67b3d25 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m @@ -28,7 +28,7 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector) B = sign(LocalMagneticField(1:3) * WaveVectorEndPoint(i,1:3)') * LocalMagneticField(4); - ZeemanShift = this.LandegFactor * this.MagneticSubLevel * Helper.PhysicsConstants.BohrMagneton * Helper.PhysicsConstants.PlanckConstantReduced * B; + ZeemanShift = this.LandegFactor * this.MagneticSubLevel * Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced * B; DopplerShift = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveVector; diff --git a/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m b/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m new file mode 100644 index 0000000..243c0b2 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m @@ -0,0 +1,52 @@ +function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ... + SecondParameterName, SecondParameterArray, varargin) + + p = inputParser; + p.KeepUnmatched = true; + + addRequired(p, 'ClassObject' , @isobject); + addRequired(p, 'FirstParameterName' , @ischar); + addRequired(p, 'FirstParameterArray' , @isvector); + addRequired(p, 'SecondParameterName' , @ischar); + addRequired(p, 'SecondParameterArray', @isvector); + + addParameter(p, 'ChangeRelatedParameter', false, @islogical); + 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); + + + p.parse(this, FirstParameterName, FirstParameterArray, ... + SecondParameterName, SecondParameterArray, varargin{:}) + + FirstParameterName = p.Results.FirstParameterName; + FirstParameterArray = p.Results.FirstParameterArray; + SecondParameterName = p.Results.SecondParameterName; + SecondParameterArray = p.Results.SecondParameterArray; + ChangeRelatedParameter = p.Results.ChangeRelatedParameter; + Order = p.Results.Order; + RelatedParameterName = p.Results.RelatedParameterName; + RelatedParameterArray = p.Results.RelatedParameterArray; + + NumberOfPointsForFirstParam = length(FirstParameterArray); + NumberOfPointsForSecondParam = length(SecondParameterArray); + + for i=1:NumberOfPointsForFirstParam + eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i))); + if ChangeRelatedParameter && Order == 1 + eval(sprintf('OptionsStruct.%s = %d;', RelatedParameterName, RelatedParameterArray(i))); + end + + for j=1:NumberOfPointsForSecondParam + eval(sprintf('OptionsStruct.%s = %d;', SecondParameterName, SecondParameterArray(j))); + if ChangeRelatedParameter && Order == 2 + eval(sprintf('OptionsStruct.%s = %d;', RelatedParameterName, RelatedParameterArray(j))); + end + options = Helper.convertstruct2cell(OptionsStruct); + this.setInitialConditions(options{:}); + tic + [LoadingRate, ~] = this.runSimulation(); + LoadingRateArray(i,j) = LoadingRate; + end + end +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m b/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m index 129f17f..4353229 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m +++ b/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m @@ -1,39 +1,39 @@ -function ret = initialVelocitySampling(this, VelocityCutoff, AngularDistribution, NormalizationConstant) +function ret = initialVelocitySampling(this) n = this.NumberOfAtoms; ret = zeros(n,3); - SampledVelocityMagnitude =zeros(n,1); - SampledPolarAngle =zeros(n,1); - SampledAzimuthalAngle =zeros(n,1); + SampledVelocityMagnitude = zeros(n,1); + SampledPolarAngle = zeros(n,1); + SampledAzimuthalAngle = zeros(n,1); MostProbableVelocity = sqrt((3 * Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperature) / Helper.PhysicsConstants.Dy164Mass); % For v * f(v) distribution - if MostProbableVelocity > VelocityCutoff - MaximumVelocityAllowed = VelocityCutoff; + if MostProbableVelocity > this.VelocityCutoff + MaximumVelocityAllowed = this.VelocityCutoff; else MaximumVelocityAllowed = MostProbableVelocity; end ProbabilityOfMaximumVelocityAllowed = this.velocityDistributionFunction(MaximumVelocityAllowed); - ProbabilityOfMaximumDivergenceAngleAllowed = NormalizationConstant * max(AngularDistribution); + ProbabilityOfMaximumDivergenceAngleAllowed = 0.98 * this.NormalizationConstantForAngularDistribution * max(this.AngularDistribution .* sin(this.ThetaArray)); parfor i = 1:n % Rejection Sampling of speed y = ProbabilityOfMaximumVelocityAllowed * rand(1); - x = VelocityCutoff * rand(1); + x = this.VelocityCutoff * rand(1); while y > this.velocityDistributionFunction(x) %As long as this loop condition is satisfied, reject the corresponding x value y = ProbabilityOfMaximumVelocityAllowed * rand(1); - x = VelocityCutoff * rand(1); + x = this.VelocityCutoff * rand(1); end SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample % Rejection Sampling of polar angle - w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); z = this.MOTExitDivergence * rand(1); + w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); - while w > (NormalizationConstant * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value - w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); + while w > (this.NormalizationConstantForAngularDistribution * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value z = this.MOTExitDivergence * rand(1); + w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); end SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample diff --git a/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m index 82bdf49..f15b8d4 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m +++ b/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m @@ -6,14 +6,13 @@ function reinitializeSimulator(this) this.NozzleRadius = 2.50e-3; this.Beta = 2 * (this.NozzleRadius/this.NozzleLength); this.ApertureCut = max(2.5e-3,this.NozzleRadius); - this.OvenDistance = ((25+12.5)*1e-3 + (this.NozzleRadius + this.ApertureCut)) / tan(15/360 * 2 * pi); + this.OvenDistance = (25+12.5)*1e-3 + (this.NozzleRadius + this.ApertureCut) / tan(15/360 * 2 * pi); % Distance between the nozzle and the 2-D MOT chamber center % 25 is the beam radius/sqrt(2) % 12.5 is the radius of the oven % 15 eg is the angle between the 2-D MOT chamber center and the nozzle this.OvenTemperature = 1000; % Temperature in Celsius this.MOTDistance = 320e-3; % Distance between the 2-D MOT the 3-D MOT - this.MagneticGradient = 0.425; % T/m this.BlueWaveVector = 2*pi/pc.BlueWavelength; this.BlueSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.BlueLinewidth/3/(pc.BlueWavelength)^3/10; this.OrangeWaveVector = 2*pi/pc.OrangeWavelength; @@ -25,8 +24,8 @@ function reinitializeSimulator(this) this.MOTExitDivergence = 0.016; % The limitation angle between 2D-MOT and 3D-MOT this.TotalPower = 0.4; this.OrangeBeamRadius = 1.2e-03; - this.PushBeamRadius = 0; - this.PushBeamDistance = 0; - this.DistanceBetweenPushBeamAnd3DMOTCenter = 1; + this.PushBeamRadius = 0.000; + this.PushBeamDistance = 0.32; + this.DistanceBetweenPushBeamAnd3DMOTCenter = 0; this.ZeemanSlowerBeamRadius = 1; 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 6c4e0f5..809220a 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m +++ b/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m @@ -32,7 +32,8 @@ function [LoadingRate, StandardError] = runSimulation(this) 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; @@ -40,7 +41,7 @@ function [LoadingRate, StandardError] = runSimulation(this) % - sampling the position distribution this.InitialPositions = this.initialPositionSampling(); % - sampling the velocity distribution - this.InitialVelocities = this.initialVelocitySampling(this.VelocityCutoff, this.AngularDistribution, this.NormalizationConstantForAngularDistribution); + this.InitialVelocities = this.initialVelocitySampling(); %% Solve ODE progressbar = Helper.parforNotifications(); @@ -52,15 +53,10 @@ function [LoadingRate, StandardError] = runSimulation(this) Velocities = this.InitialVelocities; parfor Index = 1:n ret = this.solver(Positions(Index,:), Velocities(Index,:)); - FinalDynamicalQuantities(Index,:) = ret(2); + FinalDynamicalQuantities(Index,:) = ret(end,:); progressbar.PB_iterate(); end clear Index %% Calculate the Loading Rate [LoadingRate, StandardError] = this.calculateLoadingRate(this.CaptureFraction, this.ClausingFactor, FinalDynamicalQuantities); - %% Save - %if this.DoSave - - %end - 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 64ffe39..f5edc79 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m +++ b/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m @@ -6,7 +6,7 @@ function setInitialConditions(this, varargin) @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'BluePower', 200e-3,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'BlueDetuning', -1.92857*Helper.PhysicsConstants.BlueLinewidth,... + addParameter(p, 'BlueDetuning', -1.5*Helper.PhysicsConstants.BlueLinewidth,... @(x) assert(isnumeric(x) && isscalar(x))); addParameter(p, 'BlueBeamWaist', 10e-3,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @@ -34,8 +34,10 @@ function setInitialConditions(this, varargin) @(x) assert(isnumeric(x) && isscalar(x))); addParameter(p, 'ZeemanSlowerBeamWaist', 7e-3,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - - p.parse(varargin{:}); + addParameter(p, 'MagneticGradient', 0.425,... % T/m + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + + p.parse(varargin{:}); this.NumberOfAtoms = p.Results.NumberOfAtoms; this.BluePower = p.Results.BluePower; @@ -52,7 +54,8 @@ function setInitialConditions(this, varargin) this.OrangeBeamWaist = p.Results.OrangeBeamWaist; this.ZeemanSlowerBeamPower = p.Results.ZeemanSlowerBeamPower; this.ZeemanSlowerBeamDetuning = p.Results.ZeemanSlowerBeamDetuning; - this.ZeemanSlowerBeamWaist = p.Results.ZeemanSlowerBeamWaist; + this.ZeemanSlowerBeamWaist = p.Results.ZeemanSlowerBeamWaist; + this.MagneticGradient = p.Results.MagneticGradient; %% Set general parameters according to simulation mode switch this.SimulationMode @@ -78,21 +81,61 @@ function setInitialConditions(this, varargin) %% - store in struct - this.InitialParameters = struct; - this.InitialParameters.NumberOfAtoms = this.NumberOfAtoms; - this.InitialParameters.BluePower = this.BluePower; - this.InitialParameters.BlueDetuning = this.BlueDetuning; - this.InitialParameters.BlueBeamWaist = this.BlueBeamWaist; - this.InitialParameters.SidebandPower = this.SidebandPower; - this.InitialParameters.SidebandDetuning = this.SidebandDetuning; - this.InitialParameters.SidebandBeamWaist = this.SidebandBeamWaist; - this.InitialParameters.PushBeamPower = this.PushBeamPower; - this.InitialParameters.PushBeamDetuning = this.PushBeamDetuning; - this.InitialParameters.PushBeamWaist = this.PushBeamWaist; - this.InitialParameters.OrangePower = this.OrangePower; - this.InitialParameters.OrangeDetuning = this.OrangeDetuning; - this.InitialParameters.OrangeBeamWaist = this.OrangeBeamWaist; - this.InitialParameters.ZeemanSlowerBeamPower = this.ZeemanSlowerBeamPower; - this.InitialParameters.ZeemanSlowerBeamDetuning = this.ZeemanSlowerBeamDetuning; - this.InitialParameters.ZeemanSlowerBeamBeamWaist = this.ZeemanSlowerBeamWaist; + this.SimulationParameters = struct; + this.SimulationParameters.SimulationMode = this.SimulationMode; + this.SimulationParameters.TimeStep = this.TimeStep; + this.SimulationParameters.SimulationTime = this.SimulationTime; + this.SimulationParameters.NumberOfAtoms = this.NumberOfAtoms; + this.SimulationParameters.NozzleLength = this.NozzleLength; + this.SimulationParameters.NozzleRadius = this.NozzleRadius; + this.SimulationParameters.Beta = this.Beta; + this.SimulationParameters.ApertureCut = this.ApertureCut; + this.SimulationParameters.OvenDistance = this.OvenDistance; + this.SimulationParameters.OvenTemperature = this.OvenTemperature; + this.SimulationParameters.MagneticGradient = this.MagneticGradient; + this.SimulationParameters.NozzleExitDivergence = this.NozzleExitDivergence; + this.SimulationParameters.MOTExitDivergence = this.MOTExitDivergence; + this.SimulationParameters.MOTDistance = this.MOTDistance; + this.SimulationParameters.BluePower = this.BluePower; + this.SimulationParameters.BlueDetuning = this.BlueDetuning; + this.SimulationParameters.BlueBeamRadius = this.BlueBeamRadius; + this.SimulationParameters.BlueBeamWaist = this.BlueBeamWaist; + this.SimulationParameters.BlueSaturationIntensity = this.BlueSaturationIntensity; + this.SimulationParameters.OrangePower = this.OrangePower; + this.SimulationParameters.OrangeDetuning = this.OrangeDetuning; + this.SimulationParameters.OrangeBeamRadius = this.OrangeBeamRadius; + this.SimulationParameters.OrangeBeamWaist = this.OrangeBeamWaist; + this.SimulationParameters.OrangeSaturationIntensity = this.OrangeSaturationIntensity; + this.SimulationParameters.SidebandPower = this.SidebandPower; + this.SimulationParameters.SidebandDetuning = this.SidebandDetuning; + this.SimulationParameters.SidebandBeamRadius = this.SidebandBeamRadius; + this.SimulationParameters.SidebandBeamWaist = this.SidebandBeamWaist; + this.SimulationParameters.SidebandBeamSaturationIntensity = this.SidebandBeamSaturationIntensity; + this.SimulationParameters.PushBeamPower = this.PushBeamPower; + this.SimulationParameters.PushBeamDetuning = this.PushBeamDetuning; + this.SimulationParameters.PushBeamRadius = this.PushBeamRadius; + this.SimulationParameters.PushBeamWaist = this.PushBeamWaist; + this.SimulationParameters.PushBeamDistance = this.PushBeamDistance; + this.SimulationParameters.DistanceBetweenPushBeamAnd3DMOTCenter = this.DistanceBetweenPushBeamAnd3DMOTCenter; + this.SimulationParameters.PushBeamSaturationIntensity = this.PushBeamSaturationIntensity; + this.SimulationParameters.TotalPower = this.TotalPower; + this.SimulationParameters.LandegFactor = this.LandegFactor; + this.SimulationParameters.MagneticSubLevel = this.MagneticSubLevel; + this.SimulationParameters.CoolingBeamSaturationParameter = this.CoolingBeamSaturationParameter; + this.SimulationParameters.SidebandSaturationParameter = this.SidebandSaturationParameter; + this.SimulationParameters.PushBeamSaturationParameter = this.PushBeamSaturationParameter; + this.SimulationParameters.OvenTemperatureinKelvin = this.OvenTemperatureinKelvin; + this.SimulationParameters.AverageVelocity = this.AverageVelocity; + this.SimulationParameters.AtomicBeamDensity = this.AtomicBeamDensity; + this.SimulationParameters.MeanFreePath = this.MeanFreePath; + this.SimulationParameters.CollisionTime = this.CollisionTime; + + if strcmpi(this.SimulationMode, '3D') + this.SimulationParameters.ZeemanSlowerBeamPower = this.ZeemanSlowerBeamPower; + this.SimulationParameters.ZeemanSlowerBeamDetuning = this.ZeemanSlowerBeamDetuning; + this.SimulationParameters.ZeemanSlowerBeamRadius = this.ZeemanSlowerBeamRadius; + this.SimulationParameters.ZeemanSlowerBeamBeamWaist = this.ZeemanSlowerBeamWaist; + this.SimulationParameters.ZeemanSlowerBeamSaturationIntensity = this.ZeemanSlowerBeamSaturationIntensity; + this.SimulationParameters.ZeemanSlowerBeamSaturationParameter = this.ZeemanSlowerBeamSaturationParameter; + end end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/solver.m b/MOT Capture Process Simulation/@MOTSimulator/solver.m index b0c9dfd..49e5451 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/solver.m +++ b/MOT Capture Process Simulation/@MOTSimulator/solver.m @@ -1,6 +1,6 @@ function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity) if this.Gravity - g = [0,0,- -Helper.PhysicsConstants.GravitationalAcceleration]; + g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration]; else g = 0; end @@ -9,7 +9,7 @@ function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPo collision = rand(1); CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime); - if collision >= CollisionProbability || this.CollisionTime == -500 % -500 is a flag for skipping the background collision + if collision >= CollisionProbability || this.AtomicBeamCollision == false %flag for skipping the background collision ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9); diff --git a/MOT Capture Process Simulation/test_MOTSimulator.m b/MOT Capture Process Simulation/test_MOTSimulator.m index 90a28fd..e0af22c 100644 --- a/MOT Capture Process Simulation/test_MOTSimulator.m +++ b/MOT Capture Process Simulation/test_MOTSimulator.m @@ -4,12 +4,13 @@ clc OptionsStruct = struct; OptionsStruct.SimulationMode = '2D'; -OptionsStruct.TimeStep = 50e-06; -OptionsStruct.SimulationTime = 04e-03; -OptionsStruct.SpontaneousEmission = true; -OptionsStruct.Sideband = true; +OptionsStruct.TimeStep = 50e-06; % in s +OptionsStruct.SimulationTime = 4e-03; % in s +OptionsStruct.SpontaneousEmission = false; +OptionsStruct.Sideband = false; OptionsStruct.ZeemanSlowerBeam = false; OptionsStruct.Gravity = true; +OptionsStruct.AtomicBeamCollision = true; OptionsStruct.DebugMode = false; OptionsStruct.SaveData = false; options = Helper.convertstruct2cell(OptionsStruct); @@ -24,8 +25,8 @@ Simulator.setInitialConditions(); OptionsStruct = struct; OptionsStruct.NumberOfAtoms = 5000; OptionsStruct.BluePower = 0.2; % in W -OptionsStruct.BlueDetuning = -2 * Helper.PhysicsConstants.BlueLinewidth; % in Hz -OptionsStruct.BlueBeamWaist = 0.010; % in m +OptionsStruct.BlueDetuning = -1.5 * Helper.PhysicsConstants.BlueLinewidth; % in Hz +OptionsStruct.BlueBeamWaist = 0.016; % in m OptionsStruct.SidebandPower = 0.2; OptionsStruct.SidebandDetuning = -3 * Helper.PhysicsConstants.BlueLinewidth; % in Hz OptionsStruct.SidebandBeamWaist = 0.010; % in m @@ -49,20 +50,53 @@ YAxisRange = [-5 5]; ZAxisRange = [-5 5]; Plotting.visualizeMagneticField(Simulator, XAxisRange, YAxisRange, ZAxisRange) +%% - Plot MFP & VP for different temperatures +TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius +Plotting.plotMeanFreePathAndVapourPressure(TemperatureinCelsius) + +%% - Plot the Free Molecular Flux for different temperatures +Temperature = [900, 950, 1000, 1050, 1100]; % Temperature' +Plotting.plotFreeMolecularFluxVsTemp(Simulator,Temperature) + +%% - Plot Angular Distribution for different Beta +Beta = [0.5, 0.1 , 0.05, 0.02, 0.01]; %Beta = 2 * radius / length of the tube +Plotting.plotAngularDistributionForDifferentBeta(Simulator, Beta) + +%% - Plot Capture Velocity + +Simulator.setInitialConditions(); +Simulator.SimulationTime = 15*10^(-4); +Simulator.TimeStep = 1.5*10^(-6); +Simulator.MOTExitDivergence = 15/360*2*pi; +Simulator.MagneticGradient = 0.4; + +Plotting.plotCaptureVelocityVsAngle(Simulator); + +%% - Plot Phase Space with Acceleration Field + +Simulator.NumberOfAtoms = 100; +MaximumVelocity = 150; +NumberOfBins = 200; %Along each axis +IncidentAtomDirection = 0*2*pi/360; +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 -NumberOfPointsForFirstParam = 10; %iterations of the simulation -NumberOfPointsForSecondParam = 10; +NumberOfPointsForFirstParam = 2; %iterations of the simulation +NumberOfPointsForSecondParam = 2; LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); % Scan Sideband Detuning and Power Ratio DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam); PowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam); +OptionsStruct.NumberOfAtoms = 5000; + tStart = tic; for i=1:NumberOfPointsForFirstParam OptionsStruct.SidebandDetuning = DetuningArray(i) * Helper.PhysicsConstants.BlueLinewidth; @@ -79,6 +113,30 @@ end tEnd = toc(tStart); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); +%% TWO-PARAMETER SCAN FUNCTION + +NumberOfPointsForFirstParam = 2; %iterations of the simulation +NumberOfPointsForSecondParam = 2; +Simulator.NumberOfAtoms = 5000; + +% Scan Sideband Detuning and Power Ratio +DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth; +SidebandPowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam) * Simulator.TotalPower; +BluePowerArray = Simulator.TotalPower - SidebandPowerArray; + +OptionsStruct = struct; +OptionsStruct.ChangeRelatedParameter = true; +OptionsStruct.Order = 2; %Change after first parameter = 1, Change after second parameter = 2 +OptionsStruct.RelatedParameterName = 'BluePower'; +OptionsStruct.RelatedParameterArray = BluePowerArray; +options = Helper.convertstruct2cell(OptionsStruct); + +tStart = tic; +LoadingRateArray = 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; @@ -87,9 +145,9 @@ QuantityOfInterestArray = LoadingRateArray; OptionsStruct = struct; OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1; -OptionsStruct.XLabelString = 'Detuning (\Delta/\Gamma)'; +OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)'; OptionsStruct.RescalingFactorForSecondParameter = 1000; -OptionsStruct.YLabelString = 'Sideband Beam Power (mW)'; +OptionsStruct.YLabelString = 'Sideband Power (mW)'; OptionsStruct.ZLabelString = 'Loading rate (atoms/s)'; OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);