From a9b7045c57e782b248ed76f543bac7d5c1eb80c1 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 14 Jul 2021 20:23:00 +0200 Subject: [PATCH] MAJOR UPDATE: Complete refactoring of the MOT simulation code to better reflect OOP. --- .../+Helper => +Helper}/ImageSelection.class | Bin .../+Helper => +Helper}/ImageSelection.java | 0 .../+Helper => +Helper}/PhysicsConstants.m | 14 +- .../bringFiguresWithTagInForeground.m | 0 .../calculateDistanceFromPointToLine.m | 0 .../+Helper => +Helper}/convertstruct2cell.m | 0 .../findAllZeroCrossings.m | 0 .../+Helper => +Helper}/getFigureByTag.m | 0 .../+Helper => +Helper}/ode5.m | 0 .../+Helper => +Helper}/onenoteccdata.m | 0 .../+Helper => +Helper}/parforNotifications.m | 0 .../+Helper => +Helper}/screencapture.m | 0 .../plotAngularDistributionForDifferentBeta.m | 3 - .../plotCaptureVelocityVsAngle.m | 6 +- .../plotConfidenceIntervalRegion.m | 0 +Plotter/plotDynamicalQuantities.m | 85 +++ .../plotFreeMolecularFluxVsTemp.m | 11 +- .../plotInitialVeloctiySamplingVsAngle.m | 14 +- .../plotMeanFreePathAndVapourPressureVsTemp.m | 0 .../plotPhaseSpaceWithAccelerationField.m | 19 +- .../plotPositionAndVelocitySampling.m | 0 .../plotResultForOneParameterScan.m | 0 .../plotResultForTwoParameterScan.m | 0 .../visualizeMagneticField.m | 0 +Simulator/@Beams/Beams.m | 203 ++++++ .../@MOTCaptureProcess/MOTCaptureProcess.m | 173 +++++ +Simulator/@Oven/Oven.m | 177 +++++ .../@Oven}/angularDistributionFunction.m | 0 .../@Oven}/calculateClausingFactor.m | 0 .../@Oven}/calculateFreeMolecularRegimeFlux.m | 0 .../@Oven}/calculateReducedClausingFactor.m | 2 +- +Simulator/@Oven/initialPositionSampling.m | 7 + +Simulator/@Oven/initialVelocitySampling.m | 59 ++ .../@Oven}/velocityDistributionFunction.m | 0 .../@TwoDimensionalMOT/TwoDimensionalMOT.m | 169 +++++ .../accelerationDueToPushBeam.m | 35 + ...elerationDueToSpontaneousEmissionProcess.m | 0 .../bootstrapErrorEstimation.m | 7 +- .../calculateCaptureVelocity.m | 22 + .../@TwoDimensionalMOT/calculateLoadingRate.m | 29 + .../calculateLocalSaturationIntensity.m | 0 .../calculateTotalAcceleration.m | 49 +- .../computeCollisionProbability.m | 4 +- .../computeTimeSpentInInteractionRegion.m | 4 +- .../@TwoDimensionalMOT}/exitCondition.m | 4 +- .../@TwoDimensionalMOT/magneticFieldForMOT.m | 8 + .../@TwoDimensionalMOT}/runSimulation.m | 10 +- .../@TwoDimensionalMOT}/solver.m | 1 - .../+Plotting/plotDynamicalQuantities.m | 107 --- .../@MOTSimulator/MOTSimulator.m | 654 ------------------ .../@MOTSimulator/accelerationDueToPushBeam.m | 27 - .../@MOTSimulator/calculateCaptureVelocity.m | 28 - .../@MOTSimulator/calculateLoadingRate.m | 39 -- .../@MOTSimulator/doOneParameterScan.m | 67 -- .../@MOTSimulator/doTwoParameterScan.m | 82 --- .../@MOTSimulator/initialPositionSampling.m | 12 - .../@MOTSimulator/initialVelocitySampling.m | 67 -- .../@MOTSimulator/magneticFieldForMOT.m | 14 - .../@MOTSimulator/reinitializeSimulator.m | 35 - .../@MOTSimulator/setInitialConditions.m | 137 ---- .../test_MOTSimulator.m | 188 ----- test_MOTCaptureProcessSimulation.m | 188 +++++ 62 files changed, 1231 insertions(+), 1529 deletions(-) rename {MOT Capture Process Simulation/+Helper => +Helper}/ImageSelection.class (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/ImageSelection.java (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/PhysicsConstants.m (84%) rename {MOT Capture Process Simulation/+Helper => +Helper}/bringFiguresWithTagInForeground.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/calculateDistanceFromPointToLine.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/convertstruct2cell.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/findAllZeroCrossings.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/getFigureByTag.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/ode5.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/onenoteccdata.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/parforNotifications.m (100%) rename {MOT Capture Process Simulation/+Helper => +Helper}/screencapture.m (100%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotAngularDistributionForDifferentBeta.m (96%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotCaptureVelocityVsAngle.m (79%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotConfidenceIntervalRegion.m (100%) create mode 100644 +Plotter/plotDynamicalQuantities.m rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotFreeMolecularFluxVsTemp.m (89%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotInitialVeloctiySamplingVsAngle.m (84%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotMeanFreePathAndVapourPressureVsTemp.m (100%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotPhaseSpaceWithAccelerationField.m (76%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotPositionAndVelocitySampling.m (100%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotResultForOneParameterScan.m (100%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/plotResultForTwoParameterScan.m (100%) rename {MOT Capture Process Simulation/+Plotting => +Plotter}/visualizeMagneticField.m (100%) create mode 100644 +Simulator/@Beams/Beams.m create mode 100644 +Simulator/@MOTCaptureProcess/MOTCaptureProcess.m create mode 100644 +Simulator/@Oven/Oven.m rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@Oven}/angularDistributionFunction.m (100%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@Oven}/calculateClausingFactor.m (100%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@Oven}/calculateFreeMolecularRegimeFlux.m (100%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@Oven}/calculateReducedClausingFactor.m (94%) create mode 100644 +Simulator/@Oven/initialPositionSampling.m create mode 100644 +Simulator/@Oven/initialVelocitySampling.m rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@Oven}/velocityDistributionFunction.m (100%) create mode 100644 +Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m create mode 100644 +Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/accelerationDueToSpontaneousEmissionProcess.m (100%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/bootstrapErrorEstimation.m (90%) create mode 100644 +Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m create mode 100644 +Simulator/@TwoDimensionalMOT/calculateLoadingRate.m rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/calculateLocalSaturationIntensity.m (100%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/calculateTotalAcceleration.m (58%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/computeCollisionProbability.m (53%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/computeTimeSpentInInteractionRegion.m (81%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/exitCondition.m (57%) create mode 100644 +Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/runSimulation.m (77%) rename {MOT Capture Process Simulation/@MOTSimulator => +Simulator/@TwoDimensionalMOT}/solver.m (99%) delete mode 100644 MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m delete mode 100644 MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m delete mode 100644 MOT Capture Process Simulation/test_MOTSimulator.m create mode 100644 test_MOTCaptureProcessSimulation.m diff --git a/MOT Capture Process Simulation/+Helper/ImageSelection.class b/+Helper/ImageSelection.class similarity index 100% rename from MOT Capture Process Simulation/+Helper/ImageSelection.class rename to +Helper/ImageSelection.class diff --git a/MOT Capture Process Simulation/+Helper/ImageSelection.java b/+Helper/ImageSelection.java similarity index 100% rename from MOT Capture Process Simulation/+Helper/ImageSelection.java rename to +Helper/ImageSelection.java diff --git a/MOT Capture Process Simulation/+Helper/PhysicsConstants.m b/+Helper/PhysicsConstants.m similarity index 84% rename from MOT Capture Process Simulation/+Helper/PhysicsConstants.m rename to +Helper/PhysicsConstants.m index 8f8cc49..f62dd0a 100644 --- a/MOT Capture Process Simulation/+Helper/PhysicsConstants.m +++ b/+Helper/PhysicsConstants.m @@ -28,14 +28,14 @@ classdef PhysicsConstants < handle 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; + BlueLinewidth = 1/4.94e-9; + RedWavelength = 626.086e-9; + RedLandegFactor = 1.29; + RedLifetime = 1.2e-6; + RedLinewidth = 1/1.2e-6; PushBeamWaveLength = 626.086e-9; - PushBeamLinewidth = 8.5e5; + PushBeamLifetime = 1.2e-6; + PushBeamLinewidth = 1/1.2e-6; end methods diff --git a/MOT Capture Process Simulation/+Helper/bringFiguresWithTagInForeground.m b/+Helper/bringFiguresWithTagInForeground.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/bringFiguresWithTagInForeground.m rename to +Helper/bringFiguresWithTagInForeground.m diff --git a/MOT Capture Process Simulation/+Helper/calculateDistanceFromPointToLine.m b/+Helper/calculateDistanceFromPointToLine.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/calculateDistanceFromPointToLine.m rename to +Helper/calculateDistanceFromPointToLine.m diff --git a/MOT Capture Process Simulation/+Helper/convertstruct2cell.m b/+Helper/convertstruct2cell.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/convertstruct2cell.m rename to +Helper/convertstruct2cell.m diff --git a/MOT Capture Process Simulation/+Helper/findAllZeroCrossings.m b/+Helper/findAllZeroCrossings.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/findAllZeroCrossings.m rename to +Helper/findAllZeroCrossings.m diff --git a/MOT Capture Process Simulation/+Helper/getFigureByTag.m b/+Helper/getFigureByTag.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/getFigureByTag.m rename to +Helper/getFigureByTag.m diff --git a/MOT Capture Process Simulation/+Helper/ode5.m b/+Helper/ode5.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/ode5.m rename to +Helper/ode5.m diff --git a/MOT Capture Process Simulation/+Helper/onenoteccdata.m b/+Helper/onenoteccdata.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/onenoteccdata.m rename to +Helper/onenoteccdata.m diff --git a/MOT Capture Process Simulation/+Helper/parforNotifications.m b/+Helper/parforNotifications.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/parforNotifications.m rename to +Helper/parforNotifications.m diff --git a/MOT Capture Process Simulation/+Helper/screencapture.m b/+Helper/screencapture.m similarity index 100% rename from MOT Capture Process Simulation/+Helper/screencapture.m rename to +Helper/screencapture.m diff --git a/MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m b/+Plotter/plotAngularDistributionForDifferentBeta.m similarity index 96% rename from MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m rename to +Plotter/plotAngularDistributionForDifferentBeta.m index 3e73879..3812724 100644 --- a/MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m +++ b/+Plotter/plotAngularDistributionForDifferentBeta.m @@ -15,8 +15,6 @@ function plotAngularDistributionForDifferentBeta(obj, Beta) hold on - obj.reinitializeSimulator(); - SimulationBeta = obj.Beta; if ~ismember(SimulationBeta, Beta) @@ -38,7 +36,6 @@ function plotAngularDistributionForDifferentBeta(obj, 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) diff --git a/MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m b/+Plotter/plotCaptureVelocityVsAngle.m similarity index 79% rename from MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m rename to +Plotter/plotCaptureVelocityVsAngle.m index 7352270..92c9306 100644 --- a/MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m +++ b/+Plotter/plotCaptureVelocityVsAngle.m @@ -1,10 +1,10 @@ -function plotCaptureVelocityVsAngle(obj) +function plotCaptureVelocityVsAngle(OvenObj, MOTObj) - theta = linspace(0, obj.NozzleExitDivergence, 1000); + theta = linspace(0, OvenObj.ExitDivergence, 1000); 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))]); + CaptureVelocity(i,:) = MOTObj.calculateCaptureVelocity(OvenObj, [-OvenObj.OvenDistance,0,0], [cos(theta(i)),0,sin(theta(i))]); end f_h = Helper.getFigureByTag('Capture Velocity'); diff --git a/MOT Capture Process Simulation/+Plotting/plotConfidenceIntervalRegion.m b/+Plotter/plotConfidenceIntervalRegion.m similarity index 100% rename from MOT Capture Process Simulation/+Plotting/plotConfidenceIntervalRegion.m rename to +Plotter/plotConfidenceIntervalRegion.m diff --git a/+Plotter/plotDynamicalQuantities.m b/+Plotter/plotDynamicalQuantities.m new file mode 100644 index 0000000..dad74a4 --- /dev/null +++ b/+Plotter/plotDynamicalQuantities.m @@ -0,0 +1,85 @@ +function plotDynamicalQuantities(OvenObj, MOTObj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, varargin) + + p = inputParser; + p.KeepUnmatched = true; + + addRequired(p, 'OvenObject' , @isobject); + addRequired(p, 'MOTObject' , @isobject); + addRequired(p, 'MaximumVelocity' , @(x) assert(isnumeric(x) && isscalar(x))); + addRequired(p, 'IncidentAtomDirection' , @(x) assert(isnumeric(x) && isscalar(x))); + addRequired(p, 'IncidentAtomPosition' , @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'PlotPositions' , false, @islogical); + addParameter(p, 'PlotVelocities' , false, @islogical); + + p.parse(OvenObj, MOTObj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, varargin{:}) + + MaximumVelocity = p.Results.MaximumVelocity; + IncidentAtomDirection = p.Results.IncidentAtomDirection; + IncidentAtomPosition = p.Results.IncidentAtomDirection; + PlotPositions = p.Results.PlotPositions; + PlotVelocities = p.Results.PlotVelocities; + + f_h = Helper.getFigureByTag('Trajectories 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 = 'Trajectories Plot'; + f_h.Units = 'pixels'; + + set(0,'units','pixels'); + screensize = get(0,'ScreenSize'); + f_h.Position = [[screensize(3)/50 screensize(4)/3.5] 1870 600]; + + N = MOTObj.NumberOfAtoms; + Theta = IncidentAtomDirection; + z = IncidentAtomPosition; + + L = OvenObj.OvenDistance * 2; + T = MOTObj.SimulationTime; + tau = MOTObj.TimeStep; + + Y = linspace(0,MaximumVelocity,N); + DynamicalQuantities = zeros(length(Y),int64(T/tau),6); + 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]; + DynamicalQuantities(i,:,:) = MOTObj.solver(r, v); + end + + LabelStringArrayForPositions = {'x (mm)', 'y (mm)', 'z (mm)'}; + LabelStringArrayForVelocities = {'v_x (m/s)', 'v_y (m/s)', 'v_z (m/s)'}; + colors = { [0, 0.4470, 0.7410], [0.8500, 0.3250, 0.0980], [0.4940, 0.1840, 0.5560]}; + t = linspace(0, T, T/tau) * 1e+3; + for i=1:length(Y) + for j = 1:3 + if PlotPositions + subplot(1, 3, j) + hold on + plot(t, DynamicalQuantities(i,:,j) * 1e+3,'Color', colors{j},'linewidth',1.3) + hXLabel = xlabel('Time (ms)'); + hYLabel = ylabel(LabelStringArrayForPositions{j}); + hold off + set([hXLabel, hYLabel], 'FontSize', 14); + elseif PlotVelocities + subplot(1, 3, j) + hold on + plot(t, DynamicalQuantities(i,:,j+3),'Color', colors{j},'linewidth',1.3) + hXLabel = xlabel('Time (ms)'); + hYLabel = ylabel(LabelStringArrayForVelocities{j}); + hold off + set([hXLabel, hYLabel], 'FontSize', 14); + end + end + end + + hold off + hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", MOTObj.MagneticGradient)); + set(hTitle, 'FontSize', 18); + Helper.bringFiguresWithTagInForeground(); +end + diff --git a/MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m b/+Plotter/plotFreeMolecularFluxVsTemp.m similarity index 89% rename from MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m rename to +Plotter/plotFreeMolecularFluxVsTemp.m index 1900384..46e5734 100644 --- a/MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m +++ b/+Plotter/plotFreeMolecularFluxVsTemp.m @@ -20,19 +20,18 @@ function plotFreeMolecularFluxVsTemp(obj, Temperature) obj.OvenTemperature = Temperature(i); flux = zeros(1,length(beta)); for j=1:length(beta) - obj.Beta = beta(j); + obj.NozzleLength = (2 * obj.NozzleRadius) / beta(j); [ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor(); flux(j) = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux(); end plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5) end set(gca,'yscale','log') - - obj.reinitializeSimulator(); - [ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor(); - + + obj.restoreDefaults(); + xline(obj.Beta, 'k--','Linewidth', 0.5); - fmf = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux(); + fmf = obj.ReducedClausingFactor * 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'); diff --git a/MOT Capture Process Simulation/+Plotting/plotInitialVeloctiySamplingVsAngle.m b/+Plotter/plotInitialVeloctiySamplingVsAngle.m similarity index 84% rename from MOT Capture Process Simulation/+Plotting/plotInitialVeloctiySamplingVsAngle.m rename to +Plotter/plotInitialVeloctiySamplingVsAngle.m index 11810b5..47b9479 100644 --- a/MOT Capture Process Simulation/+Plotting/plotInitialVeloctiySamplingVsAngle.m +++ b/+Plotter/plotInitialVeloctiySamplingVsAngle.m @@ -1,8 +1,6 @@ -function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins) +function plotInitialVeloctiySamplingVsAngle(obj, MOTObj, NumberOfBins) n = obj.NumberOfAtoms; - - obj.setInitialConditions(); - VelocitySamples = initialVelocitySampling(obj); + VelocitySamples = obj.initialVelocitySampling(MOTObj); Speeds = zeros(length(VelocitySamples),1); Angles = zeros(length(VelocitySamples),1); @@ -19,8 +17,8 @@ function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins) * obj.OvenTemperatureinKelvin))); c = integral(VelocityDistribution, 0, obj.VelocityCutoff); - AnglesArray = linspace(0, obj.NozzleExitDivergence,1000); - AngleBinSize = obj.NozzleExitDivergence / NumberOfBins; + AnglesArray = linspace(0, obj.ExitDivergence,1000); + AngleBinSize = obj.ExitDivergence / NumberOfBins; dtheta = AnglesArray(2)-AnglesArray(1); AngularDistribution = zeros(1,length(AnglesArray)); d = 0; @@ -57,8 +55,8 @@ function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins) histogram(Angles .* 1e+03, NumberOfBins,'FaceAlpha',0.4, 'EdgeAlpha',0.4) hold on plot(AnglesArray .* 1e+03, (n/d * AngleBinSize .* AngularDistribution),'--', 'Linewidth',1.5) - xline(obj.NozzleExitDivergence.* 1e+03, 'k--', 'Linewidth', 1.5); - text((obj.NozzleExitDivergence - 0.5 * obj.NozzleExitDivergence).* 1e+03, max(h_1.Values) - 0.45 * max(h_1.Values), 'Maximum Allowed Divergence ----------->'); + xline(obj.ExitDivergence.* 1e+03, 'k--', 'Linewidth', 1.5); + text((obj.ExitDivergence - 0.5 * obj.ExitDivergence).* 1e+03, max(h_1.Values) - 0.45 * max(h_1.Values), 'Maximum Allowed Divergence ----------->'); hXLabel_2 = xlabel('\theta (mrad)'); hYLabel_2 = ylabel('Counts'); hold off diff --git a/MOT Capture Process Simulation/+Plotting/plotMeanFreePathAndVapourPressureVsTemp.m b/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m similarity index 100% rename from MOT Capture Process Simulation/+Plotting/plotMeanFreePathAndVapourPressureVsTemp.m rename to +Plotter/plotMeanFreePathAndVapourPressureVsTemp.m diff --git a/MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m b/+Plotter/plotPhaseSpaceWithAccelerationField.m similarity index 76% rename from MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m rename to +Plotter/plotPhaseSpaceWithAccelerationField.m index 71f915e..f3c26ea 100644 --- a/MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m +++ b/+Plotter/plotPhaseSpaceWithAccelerationField.m @@ -1,5 +1,4 @@ -function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition) - +function plotPhaseSpaceWithAccelerationField(OvenObj, MOTObj, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition) f_h = Helper.getFigureByTag('Phase Space Plot'); set(groot,'CurrentFigure',f_h); a_h = get(f_h, 'CurrentAxes'); @@ -13,22 +12,22 @@ function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVeloci screensize = get(0,'ScreenSize'); f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600]; - N = obj.NumberOfAtoms; - L = obj.OvenDistance * 2; + N = MOTObj.NumberOfAtoms; + L = OvenObj.OvenDistance * 2; Theta = IncidentAtomDirection; z = IncidentAtomPosition; - T = obj.SimulationTime; - tau = obj.TimeStep; + T = MOTObj.SimulationTime; + tau = MOTObj.TimeStep; [X,Y] = meshgrid(-L/2:L/NumberOfBins:L/2,-MaximumVelocity:2*MaximumVelocity/NumberOfBins:MaximumVelocity); a=zeros(NumberOfBins+1,NumberOfBins+1,3); - obj.setInitialConditions(); + MOTObj.restoreDefaults(); 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)]); + a(i,j,:) = MOTObj.calculateTotalAcceleration([X(1,i), 0, z], [Y(j,1)*cos(Theta),0,Y(j,1)*sin(Theta)]); end end for i=1:length(X) @@ -57,7 +56,7 @@ function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVeloci vz = Y(i)*sin(Theta); r = [x,0,z]; v = [vx,0,vz]; - DynamicalQuantities(i,:,:) = obj.solver(r, v); + DynamicalQuantities(i,:,:) = MOTObj.solver(r, v); end hold on @@ -70,7 +69,7 @@ function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVeloci hXLabel = xlabel('Position: Along the x-axis (m)'); hYLabel = ylabel('Velocity (m/s)'); - hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient)); + hTitle = sgtitle(sprintf("Magnetic gradient = %.3f T/m", MOTObj.MagneticGradient)); set([hXLabel, hYLabel] , ... 'FontSize' , 14 ); set( hTitle , ... diff --git a/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m b/+Plotter/plotPositionAndVelocitySampling.m similarity index 100% rename from MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m rename to +Plotter/plotPositionAndVelocitySampling.m diff --git a/MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m b/+Plotter/plotResultForOneParameterScan.m similarity index 100% rename from MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m rename to +Plotter/plotResultForOneParameterScan.m diff --git a/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m b/+Plotter/plotResultForTwoParameterScan.m similarity index 100% rename from MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m rename to +Plotter/plotResultForTwoParameterScan.m diff --git a/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m b/+Plotter/visualizeMagneticField.m similarity index 100% rename from MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m rename to +Plotter/visualizeMagneticField.m diff --git a/+Simulator/@Beams/Beams.m b/+Simulator/@Beams/Beams.m new file mode 100644 index 0000000..d225d7a --- /dev/null +++ b/+Simulator/@Beams/Beams.m @@ -0,0 +1,203 @@ +classdef Beams < handle & matlab.mixin.Copyable + + properties (Access = private) + + BlueBeamDefault = struct('Alias', 'Blue', ... + 'Power', 400e-3, ... + 'Detuning', -1.39*Helper.PhysicsConstants.BlueLinewidth, ... + 'Radius', 16.67e-3, ... + 'Waist', 15e-3, ... + 'WaveNumber',2*pi/Helper.PhysicsConstants.BlueWavelength, ... + 'Linewidth', Helper.PhysicsConstants.BlueLinewidth, ... + 'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ... + ((Helper.PhysicsConstants.PlanckConstantReduced * ... + Helper.PhysicsConstants.SpeedOfLight * ... + Helper.PhysicsConstants.BlueLinewidth) / (Helper.PhysicsConstants.BlueWavelength)^3)); + + BlueSidebandBeamDefault = struct('Alias', 'BlueSideband', ... + 'Power', 200e-3, ... + 'Detuning', -3*Helper.PhysicsConstants.BlueLinewidth, ... + 'Radius', 16.67e-3, ... + 'Waist', 10e-3, ... + 'WaveNumber',2*pi/Helper.PhysicsConstants.BlueWavelength, ... + 'Linewidth', Helper.PhysicsConstants.BlueLinewidth, ... + 'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ... + ((Helper.PhysicsConstants.PlanckConstantReduced * ... + Helper.PhysicsConstants.SpeedOfLight * ... + Helper.PhysicsConstants.BlueLinewidth) / (Helper.PhysicsConstants.BlueWavelength)^3)); + + PushBeamDefault = struct('Alias', 'Push', ... + 'Power', 10e-3 , ... + 'Detuning', 0*Helper.PhysicsConstants.RedLinewidth, ... + 'Radius', 1.2e-03, ... + 'Waist', 20.001e-3, ... + 'WaveNumber',2*pi/Helper.PhysicsConstants.RedWavelength, ... + 'Linewidth', Helper.PhysicsConstants.RedLinewidth, ... + 'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ... + ((Helper.PhysicsConstants.PlanckConstantReduced * ... + Helper.PhysicsConstants.SpeedOfLight * ... + Helper.PhysicsConstants.RedLinewidth) / (Helper.PhysicsConstants.RedWavelength)^3)); + + RedBeamDefault = struct('Alias', 'Red', ... + 'Power', 100e-3 , ... + 'Detuning', -1*Helper.PhysicsConstants.RedLinewidth, ... + 'Radius', 1.2e-3, ... + 'Waist', 12e-3 , ... + 'WaveNumber',2*pi/Helper.PhysicsConstants.RedWavelength, ... + 'Linewidth', Helper.PhysicsConstants.RedLinewidth, ... + 'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ... + ((Helper.PhysicsConstants.PlanckConstantReduced * ... + Helper.PhysicsConstants.SpeedOfLight * ... + Helper.PhysicsConstants.RedLinewidth) / (Helper.PhysicsConstants.RedWavelength)^3)); + end + + properties + Alias + Power; + Detuning; + Radius; + Waist; + WaveNumber; + SaturationIntensity; + Linewidth; + end + + properties (Dependent) + SaturationParameter; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods + %% Class Constructor + function this = Beams(BeamName) + input = inputParser; + addRequired(input,'BeamName', @ischar); + parse(input, BeamName); + this.Alias = input.Results.BeamName; + + switch this.Alias + case this.BlueBeamDefault.Alias + this.Power = this.BlueBeamDefault.Power; + this.Detuning = this.BlueBeamDefault.Detuning; + this.Radius = this.BlueBeamDefault.Radius; + this.Waist = this.BlueBeamDefault.Waist; + this.WaveNumber = this.BlueBeamDefault.WaveNumber; + this.Linewidth = this.BlueBeamDefault.Linewidth; + this.SaturationIntensity = this.BlueBeamDefault.SaturationIntensity; + case this.BlueSidebandBeamDefault.Alias + this.Power = this.BlueSidebandBeamDefault.Power; + this.Detuning = this.BlueSidebandBeamDefault.Detuning; + this.Radius = this.BlueSidebandBeamDefault.Radius; + this.Waist = this.BlueSidebandBeamDefault.Waist; + this.WaveNumber = this.BlueSidebandBeamDefault.WaveNumber; + this.Linewidth = this.BlueSidebandBeamDefault.Linewidth; + this.SaturationIntensity = this.BlueSidebandBeamDefault.SaturationIntensity; + case this.PushBeamDefault.Alias + this.Power = this.PushBeamDefault.Power; + this.Detuning = this.PushBeamDefault.Detuning; + this.Radius = this.PushBeamDefault.Radius; + this.Waist = this.PushBeamDefault.Waist; + this.WaveNumber = this.PushBeamDefault.WaveNumber; + this.Linewidth = this.PushBeamDefault.Linewidth; + this.SaturationIntensity = this.PushBeamDefault.SaturationIntensity; + case this.RedBeamDefault.Alias + this.Power = this.RedBeamDefault.Power; + this.Detuning = this.RedBeamDefault.Detuning; + this.Radius = this.RedBeamDefault.Radius; + this.Waist = this.RedBeamDefault.Waist; + this.WaveNumber = this.RedBeamDefault.WaveNumber; + this.Linewidth = this.RedBeamDefault.Linewidth; + this.SaturationIntensity = this.RedBeamDefault.SaturationIntensity; + otherwise + error('No such beam!') + end + end + end % - lifecycle + + methods + function set.Power(this,val) + this.Power = val; + end + function ret = get.Power(this) + ret = this.Power; + end + function set.Detuning(this, val) + this.Detuning = val; + end + function ret = get.Detuning(this) + ret = this.Detuning; + end + function set.Radius(this, val) + this.Radius = val; + end + function ret = get.Radius(this) + ret = this.Radius; + end + function set.Waist(this, val) + this.Waist = val; + end + function ret = get.Waist(this) + ret = this.Waist; + end + function set.WaveNumber(this, val) + this.WaveNumber = val; + end + function ret = get.WaveNumber(this) + ret = this.WaveNumber; + end + function set.SaturationIntensity(this, val) + this.SaturationIntensity = val; + end + function ret = get.SaturationIntensity(this) + ret = this.SaturationIntensity; + end + function set.Linewidth(this, val) + this.Linewidth = val; + end + function ret = get.Linewidth(this) + ret = this.Linewidth; + end + end % - setters and getters + + methods + function ret = get.SaturationParameter(this) + ret = 0.1 * (4 * this.Power) / (pi*this.Waist^2 * this.SaturationIntensity); % two beams are reflected + end + end % - getters for dependent properties + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods(Access = protected) + function cp = copyElement(this) + % Shallow copy object + cp = copyElement@matlab.mixin.Copyable(this); + + % Forces the setter to redefine the function handles to the new copied object + + pl = properties(this); + for k = 1:length(pl) + sc = superclasses(this.(pl{k})); + if any(contains(sc,{'matlab.mixin.Copyable'})) + cp.(pl{k}) = this.(pl{k}).copy(); + end + end + end + end + + methods (Static) + + % Creates an Instance of Class, ensures singleton behaviour (that there + % can only be one Instance of this class + function singleObj = getInstance(BeamName) + % Creates an Instance of Class, ensures singleton behaviour + persistent localObj; + if isempty(localObj) || ~isvalid(localObj) + localObj = Simulator.Beams(BeamName); + end + singleObj = localObj; + end + end +end diff --git a/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m b/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m new file mode 100644 index 0000000..86d1e26 --- /dev/null +++ b/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m @@ -0,0 +1,173 @@ +classdef MOTCaptureProcess < handle & matlab.mixin.Copyable + + properties (Access = public) + SimulationMode; % MOT type + TimeStep; + SimulationTime; + NumberOfAtoms; + + %Flags + SpontaneousEmission; + Sideband; + PushBeam; + ZeemanSlowerBeam; + Gravity; + BackgroundCollision; + + DebugMode; + DoSave; + SaveDirectory; + end + + properties + Beams = {}; %Contains beam objects + end % - public + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods + %% Class Destructor (Clears object) + function this = MOTCaptureProcess(varargin) + + p = inputParser; + p.KeepUnmatched = true; + addParameter(p, 'SimulationMode', '2D',... + @(x) any(strcmpi(x,{'2D','3D', 'Full'}))); + addParameter(p, 'NumberOfAtoms', 5000,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'TimeStep', 10e-06,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'SimulationTime', 3e-03,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'SpontaneousEmission', false,... + @islogical); + addParameter(p, 'Sideband', false,... + @islogical); + addParameter(p, 'PushBeam', false,... + @islogical); + addParameter(p, 'ZeemanSlowerBeam', false,... + @islogical); + addParameter(p, 'Gravity', false,... + @islogical); + addParameter(p, 'BackgroundCollision', false,... + @islogical); + addParameter(p, 'DebugMode', false,... + @islogical); + addParameter(p, 'SaveData', false,... + @islogical); + addParameter(p, 'SaveDirectory', pwd,... + @ischar); + + p.parse(varargin{:}); + + this.SimulationMode = p.Results.SimulationMode; + + this.NumberOfAtoms = p.Results.NumberOfAtoms; + this.TimeStep = p.Results.TimeStep; + this.SimulationTime = p.Results.SimulationTime; + + this.SpontaneousEmission = p.Results.SpontaneousEmission; + this.Sideband = p.Results.Sideband; + this.PushBeam = p.Results.PushBeam; + this.ZeemanSlowerBeam = p.Results.ZeemanSlowerBeam; + this.Gravity = p.Results.Gravity; + this.BackgroundCollision = p.Results.BackgroundCollision; + + this.DebugMode = p.Results.DebugMode; + this.DoSave = p.Results.SaveData; + this.SaveDirectory = p.Results.SaveDirectory; + + switch this.SimulationMode + case "2D" + this.Beams{1} = Simulator.Beams('Blue'); + this.Beams{2} = Simulator.Beams('BlueSideband'); + this.Beams{3} = Simulator.Beams('Push'); + case "3D" + this.Beams{1} = Simulator.Beams('Red'); + % Development In progress + case "Full" + % Development In progress + end + end + end % - lifecycle + + methods + function set.TimeStep(this, val) + assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!'); + this.TimeStep = val; + end + function ret = get.TimeStep(this) + 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!'); + this.SimulationTime = val; + end + function ret = get.SimulationTime(this) + ret = this.SimulationTime; + end + function set.NumberOfAtoms(this, val) + assert(val <= 20000, 'Not time efficient to compute for atom numbers larger than 20,000!'); + this.NumberOfAtoms = val; + end + function ret = get.NumberOfAtoms(this) + ret = this.NumberOfAtoms; + end + + function set.DebugMode(this, val) + this.DebugMode = val; + end + function ret = get.DebugMode(this) + ret = this.DebugMode; + end + function set.DoSave(this, val) + this.DoSave = val; + end + 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 + + end % - setters and getters + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods(Access = protected) + function cp = copyElement(this) + % Shallow copy object + cp = copyElement@matlab.mixin.Copyable(this); + + % Forces the setter to redefine the function handles to the new copied object + + pl = properties(this); + for k = 1:length(pl) + sc = superclasses(this.(pl{k})); + if any(contains(sc,{'matlab.mixin.Copyable'})) + cp.(pl{k}) = this.(pl{k}).copy(); + end + end + end + end + + methods (Static) + + % Creates an Instance of Class, ensures singleton behaviour (that there + % can only be one Instance of this class + function singleObj = getInstance(varargin) + % Creates an Instance of Class, ensures singleton behaviour + persistent localObj; + if isempty(localObj) || ~isvalid(localObj) + localObj = Simulator.MOTCaptureProcess(varargin{:}); + end + singleObj = localObj; + end + end + +end diff --git a/+Simulator/@Oven/Oven.m b/+Simulator/@Oven/Oven.m new file mode 100644 index 0000000..c343bd1 --- /dev/null +++ b/+Simulator/@Oven/Oven.m @@ -0,0 +1,177 @@ +classdef Oven < Simulator.MOTCaptureProcess & matlab.mixin.Copyable + + properties (Access = private) + + OvenDefaults = struct('NozzleLength', 60e-3, ... + 'NozzleRadius', 2.60e-3, ... + 'OvenTemperature', 1000); + end + + properties (Access = public) + NozzleLength; + NozzleRadius; + OvenTemperature; + + VelocityCutoff; + ClausingFactor; + ReducedClausingFactor; + ReducedFlux; + end + + properties (Dependent) + ExitDivergence; + Beta; + OvenDistance; + OvenTemperatureinKelvin; + AverageVelocity; + AtomicBeamDensity; + MeanFreePath; + CollisionTime; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods + function this = Oven(varargin) + this@Simulator.MOTCaptureProcess('SimulationMode', '2D', varargin{:}); + this.NozzleLength = this.OvenDefaults.NozzleLength; + this.NozzleRadius = this.OvenDefaults.NozzleRadius; + this.OvenTemperature = this.OvenDefaults.OvenTemperature; + this.ClausingFactor = this.calculateClausingFactor(); + [this.ReducedClausingFactor, ~] = this.calculateReducedClausingFactor(); + end + + function restoreDefaults(this) + this.NozzleLength = this.OvenDefaults.NozzleLength; + this.NozzleRadius = this.OvenDefaults.NozzleRadius; + this.OvenTemperature = this.OvenDefaults.OvenTemperature; + this.ClausingFactor = this.calculateClausingFactor(); + [this.ReducedClausingFactor, ~] = this.calculateReducedClausingFactor(); + end + + end % - lifecycle + + methods + function set.NozzleLength(this,val) + this.NozzleLength = val; + end + function ret = get.NozzleLength(this) + ret = this.NozzleLength; + end + function set.NozzleRadius(this,val) + this.NozzleRadius = val; + end + function ret = get.NozzleRadius(this) + ret = this.NozzleRadius; + end + function set.OvenTemperature(this,val) + this.OvenTemperature = val; + end + function ret = get.OvenTemperature(this) + ret = this.OvenTemperature; + end + function set.VelocityCutoff(this,val) + this.VelocityCutoff = val; + end + function ret = get.VelocityCutoff(this) + ret = this.VelocityCutoff; + end + function set.ClausingFactor(this,val) + this.ClausingFactor = val; + end + function ret = get.ClausingFactor(this) + ret = this.ClausingFactor; + end + function set.ReducedClausingFactor(this,val) + this.ReducedClausingFactor = val; + end + function ret = get.ReducedClausingFactor(this) + ret = this.ReducedClausingFactor; + end + function set.ReducedFlux(this,val) + this.ReducedFlux = val; + end + function ret = get.ReducedFlux(this) + ret = this.ReducedFlux; + end + end % - setters and getters + + methods + function ret = get.Beta(this) + ret = 2 * (this.NozzleRadius/this.NozzleLength); + end + + function ret = get.OvenDistance(this) + ApertureCut = max(2.5e-3,this.NozzleRadius); + ret = (25+12.5)*1e-3 + (this.NozzleRadius + ApertureCut) / tan(15/360 * 2 * pi); + end + + function ret = get.ExitDivergence(this) + Theta_Nozzle = atan((this.NozzleRadius + 0.035/sqrt(2))/this.OvenDistance); % The angle of capture region towards the oven nozzle + Theta_Aperture = 15/360*2*pi; % The limitation angle of the second aperture in the oven + ret = min(Theta_Nozzle,Theta_Aperture); + end + + function ret = get.OvenTemperatureinKelvin(this) + ret = this.OvenTemperature + Helper.PhysicsConstants.ZeroKelvin; + end + + function ret = get.AverageVelocity(this) + %See Background collision probability section in Barbiero + ret = sqrt((8 * pi * Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass)); + end + + function ret = get.AtomicBeamDensity(this) + %See Background collision probability section in Barbiero + ret = this.calculateFreeMolecularRegimeFlux() / (this.AverageVelocity * pi * (this.Beta*this.NozzleLength/2)^2); + end + + function ret = get.MeanFreePath(this) + % Cross section = pi ( 2 * Van-der-waals radius of Dy)^2; + % Van-der-waals radius of Dy = 281e-12 + %See Expected atomic flux section and Background collision probability section in Barbiero + ret = 1/(sqrt(2) * ( pi * (2*281e-12)^2) * this.AtomicBeamDensity); + end + + function ret = get.CollisionTime(this) + ret = 3 * this.MeanFreePath/this.AverageVelocity; %See Background collision probability section in Barbiero + end + + end % - getters for dependent properties + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods(Access = protected) + function cp = copyElement(this) + % Shallow copy object + cp = copyElement@matlab.mixin.Copyable(this); + + % Forces the setter to redefine the function handles to the new copied object + + pl = properties(this); + for k = 1:length(pl) + sc = superclasses(this.(pl{k})); + if any(contains(sc,{'matlab.mixin.Copyable'})) + cp.(pl{k}) = this.(pl{k}).copy(); + end + end + end + end + + methods (Static) + + % Creates an Instance of Class, ensures singleton behaviour (that there + % can only be one Instance of this class + function singleObj = getInstance(varargin) + % Creates an Instance of Class, ensures singleton behaviour + persistent localObj; + if isempty(localObj) || ~isvalid(localObj) + localObj = Simulator.Oven(); + end + singleObj = localObj; + end + end + +end diff --git a/MOT Capture Process Simulation/@MOTSimulator/angularDistributionFunction.m b/+Simulator/@Oven/angularDistributionFunction.m similarity index 100% rename from MOT Capture Process Simulation/@MOTSimulator/angularDistributionFunction.m rename to +Simulator/@Oven/angularDistributionFunction.m diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateClausingFactor.m b/+Simulator/@Oven/calculateClausingFactor.m similarity index 100% rename from MOT Capture Process Simulation/@MOTSimulator/calculateClausingFactor.m rename to +Simulator/@Oven/calculateClausingFactor.m diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m b/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m similarity index 100% rename from MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m rename to +Simulator/@Oven/calculateFreeMolecularRegimeFlux.m diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateReducedClausingFactor.m b/+Simulator/@Oven/calculateReducedClausingFactor.m similarity index 94% rename from MOT Capture Process Simulation/@MOTSimulator/calculateReducedClausingFactor.m rename to +Simulator/@Oven/calculateReducedClausingFactor.m index e994c16..919dc91 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateReducedClausingFactor.m +++ b/+Simulator/@Oven/calculateReducedClausingFactor.m @@ -9,7 +9,7 @@ function [ReducedClausingFactor, NormalizationConstantForAngularDistribution] = ReducedClausingFactor = 0; % We have to calculate the probability of an atom coming out of the oven subject to the physical constraint parfor p = 1:length(ThetaArray) % that the angle of divergence is not more than the angle subtended at the mouth of the nozzle - if ThetaArray(p) <= this.NozzleExitDivergence + if ThetaArray(p) <= this.ExitDivergence ReducedClausingFactor = ReducedClausingFactor + (2 * pi * sin(ThetaArray(p)) * AngularDistribution(p) * (ThetaArray(2)-ThetaArray(1))); end end diff --git a/+Simulator/@Oven/initialPositionSampling.m b/+Simulator/@Oven/initialPositionSampling.m new file mode 100644 index 0000000..85af394 --- /dev/null +++ b/+Simulator/@Oven/initialPositionSampling.m @@ -0,0 +1,7 @@ +function ret = initialPositionSampling(this) + n = this.NumberOfAtoms; + phi = 2 * pi * rand(n,1); + rho = this.Beta * 0.5 * this.NozzleLength * sqrt(rand(n,1)); + ret = [-this.OvenDistance * ones(n,1), rho.*cos(phi), rho.*sin(phi)]; +end + diff --git a/+Simulator/@Oven/initialVelocitySampling.m b/+Simulator/@Oven/initialVelocitySampling.m new file mode 100644 index 0000000..176dd3b --- /dev/null +++ b/+Simulator/@Oven/initialVelocitySampling.m @@ -0,0 +1,59 @@ +function ret = initialVelocitySampling(this, MOTObj) + n = this.NumberOfAtoms; + % Calculate Calculate Capture velocity --> Introduce velocity cutoff + MOTObj.CaptureVelocity = 1.05 * MOTObj.calculateCaptureVelocity(this, [-this.OvenDistance,0,0], [1,0,0]); + this.VelocityCutoff = MOTObj.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. + + [ReducedClausingFactor, NormalizationConstantForAngularDistribution] = this.calculateReducedClausingFactor(); + this.ReducedClausingFactor = ReducedClausingFactor; + + VelocityDistribution = @(velocity) sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin))^3 ... + * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ... + * this.OvenTemperatureinKelvin))); + + c = integral(VelocityDistribution, 0, this.VelocityCutoff); + + this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux(); + + ret = zeros(n,3); + 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 > this.VelocityCutoff + MaximumVelocityAllowed = this.VelocityCutoff; + else + MaximumVelocityAllowed = MostProbableVelocity; + end + NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed); + + parfor i = 1:n + % Rejection Sampling of speed + y = rand(1); + x = this.VelocityCutoff * rand(1); + while y > ((NormalizationConstantForVelocityDistribution)^-1 * this.velocityDistributionFunction(x)) %As long as this loop condition is satisfied, reject the corresponding x value + y = 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 = rand(1); + z = this.ExitDivergence * rand(1); + + while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value + w = rand(1); + z = this.ExitDivergence * rand(1); + end + SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample + + % Sampling of azimuthal angle + SampledAzimuthalAngle(i)= 2 * pi * rand(1); + + ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ... + SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))]; +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/velocityDistributionFunction.m b/+Simulator/@Oven/velocityDistributionFunction.m similarity index 100% rename from MOT Capture Process Simulation/@MOTSimulator/velocityDistributionFunction.m rename to +Simulator/@Oven/velocityDistributionFunction.m diff --git a/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m b/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m new file mode 100644 index 0000000..41ab8a7 --- /dev/null +++ b/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m @@ -0,0 +1,169 @@ +classdef TwoDimensionalMOT < Simulator.MOTCaptureProcess & matlab.mixin.Copyable + + properties (Access = private) + MagneticGradienDefault = 0.425; % T/m + ExitDivergenceDefault = 16e-3; + DistanceBetweenPushBeamAnd3DMOTCenterDefault = 0; + PushBeamDistanceDefault = 0.32; + end + + properties (Access = public) + TotalPower; + LandegFactor; + MagneticSubLevel; + MagneticGradient; + CaptureVelocity; + ExitDivergence; + DistanceBetweenPushBeamAnd3DMOTCenter; + PushBeamDistance; + TimeSpentInInteractionRegion; + ParticleDynamicalQuantities; + InitialParameters; + end + + methods + function this = TwoDimensionalMOT(varargin) + this@Simulator.MOTCaptureProcess('SimulationMode', '2D', varargin{:}); + p = inputParser; + p.KeepUnmatched = true; + addParameter(p, 'TotalPower', 0.8,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'LandegFactor', 1,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'MagneticSubLevel', 1,... + @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'MagneticGradient', 0.38,... + @(x) assert(isnumeric(x) && isscalar(x))); + + p.parse(varargin{:}); + + this.TotalPower = p.Results.TotalPower; + this.LandegFactor = p.Results.LandegFactor; + this.MagneticSubLevel = p.Results.MagneticSubLevel; + this.MagneticGradient = p.Results.MagneticGradient; + this.ExitDivergence = this.ExitDivergenceDefault; + this.DistanceBetweenPushBeamAnd3DMOTCenter = this.DistanceBetweenPushBeamAnd3DMOTCenterDefault; + this.PushBeamDistance = this.PushBeamDistanceDefault; + this.InitialParameters = struct; + this.InitialParameters.TotalPower = this.TotalPower; + this.InitialParameters.LandegFactor = this.LandegFactor; + this.InitialParameters.MagneticSubLevel= this.MagneticSubLevel; + this.InitialParameters.MagneticGradient= this.MagneticGradient; + + end + + function restoreDefaults(this) + this.TotalPower = this.InitialParameters.TotalPower; + this.LandegFactor = this.InitialParameters.LandegFactor; + this.MagneticSubLevel = this.InitialParameters.MagneticSubLevel; + this.MagneticGradient = this.InitialParameters.MagneticGradient; + this.ExitDivergence = this.ExitDivergenceDefault; + this.DistanceBetweenPushBeamAnd3DMOTCenter = this.DistanceBetweenPushBeamAnd3DMOTCenterDefault; + this.PushBeamDistance = this.PushBeamDistanceDefault; + end + end % - lifecycle + + methods + function set.TotalPower(this,val) + this.TotalPower = val; + end + function ret = get.TotalPower(this) + ret = this.TotalPower; + end + function set.LandegFactor(this,val) + this.LandegFactor = val; + end + function ret = get.LandegFactor(this) + ret = this.LandegFactor; + end + function set.MagneticSubLevel(this,val) + this.MagneticSubLevel = val; + end + function ret = get.MagneticSubLevel(this) + ret = this.MagneticSubLevel; + end + function set.MagneticGradient(this,val) + this.MagneticGradient = val; + end + function ret = get.MagneticGradient(this) + ret = this.MagneticGradient; + end + function set.CaptureVelocity(this,val) + this.CaptureVelocity = val; + end + function ret = get.CaptureVelocity(this) + ret = this.CaptureVelocity; + end + function set.ExitDivergence(this,val) + this.ExitDivergence = val; + end + function ret = get.ExitDivergence(this) + ret = this.ExitDivergence; + end + function set.DistanceBetweenPushBeamAnd3DMOTCenter(this,val) + this.DistanceBetweenPushBeamAnd3DMOTCenter = val; + end + function ret = get.DistanceBetweenPushBeamAnd3DMOTCenter(this) + ret = this.DistanceBetweenPushBeamAnd3DMOTCenter; + end + function set.PushBeamDistance(this,val) + this.PushBeamDistance = val; + end + function ret = get.PushBeamDistance(this) + ret = this.PushBeamDistance; + end + function set.TimeSpentInInteractionRegion(this,val) + this.TimeSpentInInteractionRegion = val; + end + function ret = get.TimeSpentInInteractionRegion(this) + ret = this.TimeSpentInInteractionRegion; + end + function set.ParticleDynamicalQuantities(this,val) + this.ParticleDynamicalQuantities = val; + end + function ret = get.ParticleDynamicalQuantities(this) + ret = this.ParticleDynamicalQuantities; + end + function set.InitialParameters(this,val) + this.InitialParameters = val; + end + function ret = get.InitialParameters(this) + ret = this.InitialParameters; + end + end % - setters and getters + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods(Access = protected) + function cp = copyElement(this) + % Shallow copy object + cp = copyElement@matlab.mixin.Copyable(this); + + % Forces the setter to redefine the function handles to the new copied object + + pl = properties(this); + for k = 1:length(pl) + sc = superclasses(this.(pl{k})); + if any(contains(sc,{'matlab.mixin.Copyable'})) + cp.(pl{k}) = this.(pl{k}).copy(); + end + end + end + end + + methods (Static) + + % Creates an Instance of Class, ensures singleton behaviour (that there + % can only be one Instance of this class + function singleObj = getInstance(varargin) + % Creates an Instance of Class, ensures singleton behaviour + persistent localObj; + if isempty(localObj) || ~isvalid(localObj) + localObj = Simulator.TwoDimensionalMOT(varargin{:}); + end + singleObj = localObj; + end + end + +end diff --git a/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m b/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m new file mode 100644 index 0000000..44d6485 --- /dev/null +++ b/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m @@ -0,0 +1,35 @@ +function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector) + % is the distance between the chamber center and the cross point of push beam and z-axis (along the gravity) + WaveVectorEndPoint = [0, 1, this.DistanceBetweenPushBeamAnd3DMOTCenter/this.PushBeamDistance]; + WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint); + Origin=[0,0,0]; + + PushBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), this.Beams)}; + PushBeamDetuning = PushBeamObj.Detuning; + PushBeamRadius = PushBeamObj.Radius; + PushBeamWaist = PushBeamObj.Waist; + PushBeamWaveNumber = PushBeamObj.WaveNumber; + PushBeamLinewidth = PushBeamObj.Linewidth; + PushBeamSaturationParameter = 0.25 * PushBeamObj.SaturationParameter; + + SaturationIntensity = this.calculateLocalSaturationIntensity(PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, PushBeamRadius, PushBeamWaist); + + DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * PushBeamWaveNumber; + + Detuning = PushBeamDetuning - DopplerShift; + + s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./PushBeamLinewidth).^2)); + a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(PushBeamLinewidth * 0.5); + a_push = a_sat .* (s_push/(1+s_push)); + + if this.SpontaneousEmission + a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, PushBeamLinewidth, PushBeamWaveNumber); + else + a_scatter = [0,0,0]; + end + + a_total = a_push + a_scatter; + + ret = a_total(1:3); +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m b/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m similarity index 100% rename from MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m rename to +Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m diff --git a/MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m b/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m similarity index 90% rename from MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m rename to +Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m index 731f8b4..d47cb5d 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m +++ b/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m @@ -1,5 +1,4 @@ -function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstimation(this, NumberOfLoadedAtoms) - +function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstimation(this, ovenObj, NumberOfLoadedAtoms) n = this.NumberOfAtoms; NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep); @@ -16,14 +15,14 @@ function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstima MeanLoadingRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means end - LoadingRate = mean(MeanLoadingRatioInEachSample) * this.ReducedFlux; + LoadingRate = mean(MeanLoadingRatioInEachSample) * ovenObj.ReducedFlux; Variance = 0; % Bootstrap Estimate of Variance for SampleNumber = 1:NumberOfBootsrapSamples Variance = Variance + (MeanLoadingRatioInEachSample(SampleNumber) - mean(MeanLoadingRatioInEachSample))^2; end - StandardError = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * this.ReducedFlux; + StandardError = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * ovenObj.ReducedFlux; ts = tinv([0.025 0.975],NumberOfBootsrapSamples-1); % T-Score ConfidenceInterval = LoadingRate + ts*StandardError; % 95% Confidence Intervals diff --git a/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m b/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m new file mode 100644 index 0000000..074fad7 --- /dev/null +++ b/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m @@ -0,0 +1,22 @@ +function ret = calculateCaptureVelocity(this, ovenObj, PositionVector, VelocityVector) + VelocityUnitVector = VelocityVector./norm(VelocityVector); + UpperLimit = 500; + LowerLimit = 0; + + for Index = 1:500 + InitialVelocity = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; + ParticleDynamicalQuantities = this.solver(PositionVector, InitialVelocity); + FinalPositionVector = ParticleDynamicalQuantities(end, 1:3); + if rssq(FinalPositionVector) <= ovenObj.OvenDistance + LowerLimit = 0.5 * (UpperLimit + LowerLimit); + else + UpperLimit = 0.5 * (UpperLimit + LowerLimit); + end + + if UpperLimit - LowerLimit < 1 + ret = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; + break; + end + end +end + diff --git a/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m b/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m new file mode 100644 index 0000000..2f93a97 --- /dev/null +++ b/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m @@ -0,0 +1,29 @@ +function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ovenObj) + n = this.NumberOfAtoms; + DynamicalQuantities = this.ParticleDynamicalQuantities; + NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep); + NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps); + CollisionEvents = zeros(1, n); + + % Include the stochastic process of background collisions + for AtomIndex = 1:n + this.TimeSpentInInteractionRegion(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(DynamicalQuantities(AtomIndex,:,1:3))); + CollisionEvents(AtomIndex) = this.computeCollisionProbability(ovenObj, this.TimeSpentInInteractionRegion(AtomIndex)); + end + + % Count the number of loaded atoms subject to conditions + for TimeIndex = 1:NumberOfTimeSteps + if TimeIndex ~= 1 + NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1); + end + for AtomIndex = 1:n + Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))'; + Velocity = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 4:6))'; + if this.exitCondition(ovenObj, Position, Velocity, CollisionEvents(AtomIndex)) + NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1; + end + end + end + + [LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(ovenObj, NumberOfLoadedAtoms); +end diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m b/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m similarity index 100% rename from MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m rename to +Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m b/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m similarity index 58% rename from MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m rename to +Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m index 03d343b..ee11e46 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m +++ b/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m @@ -1,5 +1,18 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector) + CoolingBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), this.Beams)}; + CoolingBeamDetuning = CoolingBeamObj.Detuning; + CoolingBeamRadius = CoolingBeamObj.Radius; + CoolingBeamWaist = CoolingBeamObj.Waist; + CoolingBeamWaveNumber = CoolingBeamObj.WaveNumber; + CoolingBeamLinewidth = CoolingBeamObj.Linewidth; + CoolingBeamSaturationParameter = CoolingBeamObj.SaturationParameter; + SidebandBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), this.Beams)}; + SidebandDetuning = SidebandBeamObj.Detuning; + SidebandBeamRadius = SidebandBeamObj.Radius; + SidebandBeamWaist = SidebandBeamObj.Waist; + SidebandSaturationParameter = SidebandBeamObj.SaturationParameter; + WaveVectorEndPoint = zeros(2,3); WaveVectorEndPoint(1,:) = [1,0,1]; WaveVectorEndPoint(1,:) = WaveVectorEndPoint(1,1:3)/norm(WaveVectorEndPoint(1,:)); @@ -10,11 +23,11 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector) Origin = [0,0,0]; % Calculate the Saturation Intensity at the specified point along its Gaussian Profile - CoolingBeamLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(0.25 * this.CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), this.CoolingBeamRadius, this.CoolingBeamWaist), ... - this.calculateLocalSaturationIntensity(0.25 * this.CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), this.CoolingBeamRadius, this.CoolingBeamWaist)]; + CoolingBeamLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(0.25 * CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), CoolingBeamRadius, CoolingBeamWaist), ... + this.calculateLocalSaturationIntensity(0.25 * CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), CoolingBeamRadius, CoolingBeamWaist)]; - SidebandLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(0.25 * this.SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), this.SidebandBeamRadius, this.SidebandBeamWaist), ... - this.calculateLocalSaturationIntensity(0.25 * this.SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), this.SidebandBeamRadius, this.SidebandBeamWaist)]; + SidebandLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(0.25 * SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), SidebandBeamRadius, SidebandBeamWaist), ... + this.calculateLocalSaturationIntensity(0.25 * SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), SidebandBeamRadius, SidebandBeamWaist)]; TotalAcceleration = zeros(1,3); @@ -29,25 +42,25 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector) ZeemanShift = this.LandegFactor * this.MagneticSubLevel * (Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced) * B; - DopplerShift = dot(WaveVectorEndPoint(i,:), VelocityVector) * this.CoolingBeamWaveNumber; + DopplerShift = dot(WaveVectorEndPoint(i,:), VelocityVector) * CoolingBeamWaveNumber; - Delta_Cooling(i*2-1) = this.CoolingBeamDetuning + DopplerShift + (ZeemanShift * Sigma(i)); - Delta_Cooling(i*2) = this.CoolingBeamDetuning - DopplerShift - (ZeemanShift * Sigma(i)); + Delta_Cooling(i*2-1) = CoolingBeamDetuning + DopplerShift + (ZeemanShift * Sigma(i)); + Delta_Cooling(i*2) = CoolingBeamDetuning - DopplerShift - (ZeemanShift * Sigma(i)); if this.Sideband - Delta_Sideband(i*2-1) = this.SidebandDetuning + DopplerShift + (ZeemanShift * Sigma(i)); - Delta_Sideband(i*2) = this.SidebandDetuning - DopplerShift - (ZeemanShift * Sigma(i)); + Delta_Sideband(i*2-1) = SidebandDetuning + DopplerShift + (ZeemanShift * Sigma(i)); + Delta_Sideband(i*2) = SidebandDetuning - DopplerShift - (ZeemanShift * Sigma(i)); end end SaturationParameter = [0,0,0,0,0,0,0,0]; for i = 1:2 - SaturationParameter(2*i-1) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2); - SaturationParameter(2*i) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i) /this.CoolingBeamLinewidth)^2); + SaturationParameter(2*i-1) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/CoolingBeamLinewidth)^2); + SaturationParameter(2*i) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i) /CoolingBeamLinewidth)^2); if this.Sideband - SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2); - SaturationParameter(2*i+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2); + SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/CoolingBeamLinewidth)^2); + SaturationParameter(2*i+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i)/CoolingBeamLinewidth)^2); end end @@ -55,13 +68,13 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector) for i = 1:2 - a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5); + a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(CoolingBeamLinewidth * 0.5); a_1 = a_sat .* (SaturationParameter(2*i-1)/(1 + TotalSaturationParameter)); a_2 = a_sat .* (SaturationParameter(2*i) /(1 + TotalSaturationParameter)); if this.SpontaneousEmission - a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ... - this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber); + a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber) + ... + this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber); else a_scattering = [0,0,0]; end @@ -72,8 +85,8 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector) if this.SpontaneousEmission a_scattering = a_scattering + ... - this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ... - this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4), TotalSaturationParameter, this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber); + this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber) + ... + this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber); else a_scattering = [0,0,0]; end diff --git a/MOT Capture Process Simulation/@MOTSimulator/computeCollisionProbability.m b/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m similarity index 53% rename from MOT Capture Process Simulation/@MOTSimulator/computeCollisionProbability.m rename to +Simulator/@TwoDimensionalMOT/computeCollisionProbability.m index 7b34fe7..39b0014 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/computeCollisionProbability.m +++ b/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m @@ -1,6 +1,6 @@ -function ret = computeCollisionProbability(this) +function ret = computeCollisionProbability(this, ovenObj, tau_2D) collision = rand(1); - CollisionProbability = 1 - exp(-this.TimeSpentInInteractionRegion/this.CollisionTime); + CollisionProbability = 1 - exp(-tau_2D/ovenObj.CollisionTime); if this.BackgroundCollision && collision <= CollisionProbability ret = true; else diff --git a/MOT Capture Process Simulation/@MOTSimulator/computeTimeSpentInInteractionRegion.m b/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m similarity index 81% rename from MOT Capture Process Simulation/@MOTSimulator/computeTimeSpentInInteractionRegion.m rename to +Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m index 4fae3f8..6e1aef0 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/computeTimeSpentInInteractionRegion.m +++ b/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m @@ -5,12 +5,12 @@ function T = computeTimeSpentInInteractionRegion(this, r) % T : gives the distribution of time spent in the interaction region % USAGE: % T = this.computeTimeSpentInInteractionRegion(r) - T = 0; + CoolingBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), this.Beams)}; NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep); for n = 1:(NumberOfTimeSteps - 1) dr = Helper.calculateDistanceFromPointToLine(r(n+1, :), [0 0 0], [0 0 1]); - if dr < this.CoolingBeamRadius + if dr < CoolingBeamObj.Radius A = 1; else A = 0; diff --git a/MOT Capture Process Simulation/@MOTSimulator/exitCondition.m b/+Simulator/@TwoDimensionalMOT/exitCondition.m similarity index 57% rename from MOT Capture Process Simulation/@MOTSimulator/exitCondition.m rename to +Simulator/@TwoDimensionalMOT/exitCondition.m index de84d2f..ac24d29 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/exitCondition.m +++ b/+Simulator/@TwoDimensionalMOT/exitCondition.m @@ -1,9 +1,9 @@ -function ret = exitCondition(this, PositionVector, VelocityVector, CollisionEvent) +function ret = exitCondition(this, ovenObj, PositionVector, VelocityVector, CollisionEvent) d = Helper.calculateDistanceFromPointToLine(PositionVector, [0 0 0], [0 1 0]); y = PositionVector(2); DivergenceAngle = (d/abs(y)); %DivergenceAngle = atan(norm(cross(PositionVector, )) / norm(dot(PositionVector, [0 1 0]))); - if (y >= 0) && (DivergenceAngle <= this.MOTExitDivergence) && (abs(VelocityVector(2))<=this.VelocityCutoff) && ~CollisionEvent + if (y >= 0) && (DivergenceAngle <= this.ExitDivergence) && (abs(VelocityVector(2))<=ovenObj.VelocityCutoff) && ~CollisionEvent ret = true; else ret = false; diff --git a/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m b/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m new file mode 100644 index 0000000..4fc4e1c --- /dev/null +++ b/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m @@ -0,0 +1,8 @@ +function ret = magneticFieldForMOT(this, r) + ret = zeros(1,4); + alpha = this.MagneticGradient; + ret(1) = r(3)*alpha; + ret(2) = 0; + ret(3) = r(1)*alpha; + ret(4) = sqrt(ret(1)^2+ret(2)^2+ret(3)^2); +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m b/+Simulator/@TwoDimensionalMOT/runSimulation.m similarity index 77% rename from MOT Capture Process Simulation/@MOTSimulator/runSimulation.m rename to +Simulator/@TwoDimensionalMOT/runSimulation.m index 3ff2a4d..36fe7c4 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m +++ b/+Simulator/@TwoDimensionalMOT/runSimulation.m @@ -1,14 +1,14 @@ -function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this) +function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this, ovenObj) %% - Sampling for initial positions and velocities % - sampling the position distribution - Positions = this.initialPositionSampling(); + Positions = ovenObj.initialPositionSampling(); % - sampling the velocity distribution - Velocities = this.initialVelocitySampling(); + Velocities = ovenObj.initialVelocitySampling(this); %% Solve ODE progressbar = Helper.parforNotifications(); - progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']); + progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating 2-D MOT capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']); % calculate the final position of the atoms DynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6); @@ -20,5 +20,5 @@ function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this) this.ParticleDynamicalQuantities = DynamicalQuantities; %% Calculate the Loading Rate - [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(); + [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ovenObj); end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/solver.m b/+Simulator/@TwoDimensionalMOT/solver.m similarity index 99% rename from MOT Capture Process Simulation/@MOTSimulator/solver.m rename to +Simulator/@TwoDimensionalMOT/solver.m index d49046e..a2c3cf1 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/solver.m +++ b/+Simulator/@TwoDimensionalMOT/solver.m @@ -1,5 +1,4 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelocity) - if this.Gravity g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration]; else diff --git a/MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m b/MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m deleted file mode 100644 index 83f22e0..0000000 --- a/MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m +++ /dev/null @@ -1,107 +0,0 @@ -function plotDynamicalQuantities(obj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition) - - f_h = Helper.getFigureByTag('Trajectories 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 = 'Trajectories Plot'; - f_h.Units = 'pixels'; - - set(0,'units','pixels'); - screensize = get(0,'ScreenSize'); - f_h.Position = [[screensize(3)/50 screensize(4)/10] 1870 850]; - - N = obj.NumberOfAtoms; - Theta = IncidentAtomDirection; - z = IncidentAtomPosition; - - obj.setInitialConditions(); - L = obj.OvenDistance * 2; - T = obj.SimulationTime; - tau = obj.TimeStep; - - Y = linspace(0,MaximumVelocity,N); - DynamicalQuantities = zeros(length(Y),int64(T/tau),6); - 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]; - DynamicalQuantities(i,:,:) = obj.solver(r, v); - end - - t = linspace(0, T, T/tau) * 1e+3; - for i=1:length(Y) - subplot(3, 3, 1) - hold on - plot(t, DynamicalQuantities(i,:,1) * 1e+3,'r','linewidth',1.3) - hXLabel_1 = xlabel('Time (ms)'); - hYLabel_1 = ylabel('x (mm)'); - hold off - subplot(3, 3, 2) - hold on - plot(t, DynamicalQuantities(i,:,2) * 1e+3,'g','linewidth',1.3) - hXLabel_2 = xlabel('Time (ms)'); - hYLabel_2 = ylabel('y (mm)'); - hold off - subplot(3, 3, 3) - hold on - plot(t, DynamicalQuantities(i,:,3) * 1e+3,'b','linewidth',1.3) - hXLabel_3 = xlabel('Time (ms)'); - hYLabel_3 = ylabel('z (m/s)'); - hold off - subplot(3, 3, 4) - hold on - plot(t, DynamicalQuantities(i,:,4),'r','linewidth',1.3) - hXLabel_4 = xlabel('Time (ms)'); - hYLabel_4 = ylabel('v_x (m/s)'); - hold off - subplot(3, 3, 5) - hold on - plot(t, DynamicalQuantities(i,:,5),'g','linewidth',1.3) - hXLabel_5 = xlabel('Time (ms)'); - hYLabel_5 = ylabel('v_y (m/s)'); - hold off - subplot(3, 3, 6) - hold on - plot(t, DynamicalQuantities(i,:,6),'b','linewidth',1.3) - hXLabel_6 = xlabel('Time (ms)'); - hYLabel_6 = ylabel('v_z (m/s)'); - hold off - subplot(3, 3, 7) - hold on - plot(DynamicalQuantities(i,:,1), DynamicalQuantities(i,:,4),'r','linewidth',1.3) - hXLabel_7 = xlabel('x (mm)'); - hYLabel_7 = ylabel('v_x (m/s)'); - xlim([-L/2 L/2]) - hold off - subplot(3, 3, 8) - hold on - plot(DynamicalQuantities(i,:,2), DynamicalQuantities(i,:,5),'g','linewidth',1.3) - hXLabel_8 = xlabel('y (mm)'); - hYLabel_8 = ylabel('v_y (m/s)'); - xlim([-1 1] * 1e-04) - hold off - subplot(3, 3, 9) - hold on - plot(DynamicalQuantities(i,:,3), DynamicalQuantities(i,:,6),'b','linewidth',1.3) - hXLabel_9 = xlabel('z (mm)'); - hYLabel_9 = ylabel('v_z (m/s)'); - xlim([-1 1] * 1e-04) - hold off - end - - hold off - hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient)); - set([hXLabel_1, hXLabel_2, hXLabel_3, hXLabel_4, hXLabel_5, hXLabel_6, hXLabel_7, hXLabel_8, hXLabel_9,... - hYLabel_1, hYLabel_2, hYLabel_3, hYLabel_4, hYLabel_5, hYLabel_6, hYLabel_7, hYLabel_8, hYLabel_9], ... - 'FontSize' , 14 ); - set( hTitle , ... - 'FontSize' , 18 ); - - Helper.bringFiguresWithTagInForeground(); -end - diff --git a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m deleted file mode 100644 index ed27b99..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m +++ /dev/null @@ -1,654 +0,0 @@ -classdef MOTSimulator < handle & matlab.mixin.Copyable - - properties (Access = public) - SimulationMode; % MOT type - TimeStep; - SimulationTime; - NumberOfAtoms; - - ParticleDynamicalQuantities; - - NozzleLength; - NozzleRadius; - Beta; - ApertureCut; - OvenDistance; - OvenTemperature; - MagneticGradient; - NozzleExitDivergence; - MOTExitDivergence; - MOTDistance; - - BluePower; - BlueDetuning; - BlueBeamRadius; - BlueBeamWaist; - BlueWaveNumber; - BlueSaturationIntensity; - - OrangePower; - OrangeDetuning; - OrangeBeamRadius; - OrangeBeamWaist; - OrangeWaveNumber; - OrangeSaturationIntensity; - - CoolingBeamPower; - CoolingBeamWaveNumber; - CoolingBeamLinewidth; - CoolingBeamDetuning; - CoolingBeamRadius; - CoolingBeamWaist; - CoolingBeamSaturationIntensity; - - SidebandPower; - SidebandDetuning; - SidebandBeamRadius; - SidebandBeamWaist; - SidebandBeamSaturationIntensity; - - PushBeamPower; - PushBeamWaveNumber; - PushBeamLinewidth; - PushBeamDetuning; - PushBeamRadius; - PushBeamWaist; - PushBeamDistance; - DistanceBetweenPushBeamAnd3DMOTCenter; - PushBeamSaturationIntensity; - - ZeemanSlowerBeamPower; - ZeemanSlowerBeamDetuning; - ZeemanSlowerBeamRadius; - ZeemanSlowerBeamWaist; - ZeemanSlowerBeamSaturationIntensity; - - TotalPower; - LandegFactor; - MagneticSubLevel; - - CaptureVelocity; - VelocityCutoff; - ClausingFactor; - ReducedClausingFactor; - ReducedFlux; - TimeSpentInInteractionRegion; - - %Flags - SpontaneousEmission; - Sideband; - PushBeam; - ZeemanSlowerBeam; - Gravity; - BackgroundCollision; - - DebugMode; - DoSave; - SaveDirectory; - - Results; - end - - properties (SetAccess = private, GetAccess = public) - SimulationParameters - end - - properties (Dependent, SetAccess = private) - CoolingBeamSaturationParameter; - SidebandSaturationParameter; - PushBeamSaturationParameter; - ZeemanSlowerBeamSaturationParameter; - OvenTemperatureinKelvin; - AverageVelocity; - AtomicBeamDensity; - MeanFreePath; - CollisionTime; - end - - methods - function s = MOTSimulator(varargin) - - p = inputParser; - p.KeepUnmatched = true; - addParameter(p, 'SimulationMode', '2D',... - @(x) any(strcmpi(x,{'2D','3D'}))); - addParameter(p, 'TimeStep', 10e-06,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'SimulationTime', 3e-03,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'SpontaneousEmission', false,... - @islogical); - addParameter(p, 'Sideband', false,... - @islogical); - addParameter(p, 'PushBeam', false,... - @islogical); - addParameter(p, 'ZeemanSlowerBeam', false,... - @islogical); - addParameter(p, 'Gravity', false,... - @islogical); - addParameter(p, 'BackgroundCollision', false,... - @islogical); - addParameter(p, 'DebugMode', false,... - @islogical); - addParameter(p, 'SaveData', false,... - @islogical); - addParameter(p, 'SaveDirectory', pwd,... - @ischar); - - p.parse(varargin{:}); - - s.SimulationMode = p.Results.SimulationMode; - - s.TimeStep = p.Results.TimeStep; - s.SimulationTime = p.Results.SimulationTime; - - s.SpontaneousEmission = p.Results.SpontaneousEmission; - s.Sideband = p.Results.Sideband; - s.PushBeam = p.Results.PushBeam; - s.ZeemanSlowerBeam = p.Results.ZeemanSlowerBeam; - s.Gravity = p.Results.Gravity; - s.BackgroundCollision = p.Results.BackgroundCollision; - - s.DebugMode = p.Results.DebugMode; - s.DoSave = p.Results.SaveData; - s.SaveDirectory = p.Results.SaveDirectory; - - - s.reinitializeSimulator(); - - poolobj = gcp('nocreate'); % Check if pool is open - if isempty(poolobj) - parpool; - end - - end - end % - lifecycle - - methods - function set.TimeStep(this, val) - assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!'); - this.TimeStep = val; - end - function ret = get.TimeStep(this) - 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!'); - this.SimulationTime = val; - end - function ret = get.SimulationTime(this) - ret = this.SimulationTime; - end - function set.NumberOfAtoms(this, val) - assert(val <= 20000, 'Not time efficient to compute for atom numbers larger than 20,000!'); - this.NumberOfAtoms = val; - end - function ret = get.NumberOfAtoms(this) - ret = this.NumberOfAtoms; - end - - function set.ParticleDynamicalQuantities(this,val) - this.ParticleDynamicalQuantities = val; - end - function ret = get.ParticleDynamicalQuantities(this) - ret = this.ParticleDynamicalQuantities; - end - - function set.NozzleLength(this,val) - this.NozzleLength = val; - end - function ret = get.NozzleLength(this) - ret = this.NozzleLength; - end - function set.NozzleRadius(this,val) - this.NozzleRadius = val; - end - function ret = get.NozzleRadius(this) - ret = this.NozzleRadius; - end - function set.Beta(this,val) - this.Beta = val; - end - function ret = get.Beta(this) - ret = this.Beta; - end - function set.ApertureCut(this,val) - this.ApertureCut = val; - end - function ret = get.ApertureCut(this) - ret = this.ApertureCut; - end - function set.OvenDistance(this,val) - this.OvenDistance = val; - end - function ret = get.OvenDistance(this) - ret = this.OvenDistance; - end - function set.OvenTemperature(this,val) - this.OvenTemperature = val; - end - function ret = get.OvenTemperature(this) - ret = this.OvenTemperature; - end - function set.MagneticGradient(this,val) - this.MagneticGradient = val; - end - function ret = get.MagneticGradient(this) - ret = this.MagneticGradient; - end - function set.NozzleExitDivergence(this,val) - this.NozzleExitDivergence = val; - end - function ret = get.NozzleExitDivergence(this) - ret = this.NozzleExitDivergence; - end - function set.MOTExitDivergence(this,val) - this.MOTExitDivergence = val; - end - function ret = get.MOTExitDivergence(this) - ret = this.MOTExitDivergence; - end - function set.MOTDistance(this,val) - this.MOTDistance = val; - end - function ret = get.MOTDistance(this) - ret = this.MOTDistance; - end - - function set.BluePower(this,val) - this.BluePower = val; - end - function ret = get.BluePower(this) - ret = this.BluePower; - end - function set.BlueDetuning(this, val) - this.BlueDetuning = val; - end - function ret = get.BlueDetuning(this) - ret = this.BlueDetuning; - end - function set.BlueBeamRadius(this, val) - this.BlueBeamRadius = val; - end - function ret = get.BlueBeamRadius(this) - ret = this.BlueBeamRadius; - end - function set.BlueBeamWaist(this, val) - this.BlueBeamWaist = val; - end - function ret = get.BlueBeamWaist(this) - ret = this.BlueBeamWaist; - end - function set.BlueWaveNumber(this, val) - this.BlueWaveNumber = val; - end - function ret = get.BlueWaveNumber(this) - ret = this.BlueWaveNumber; - end - function set.BlueSaturationIntensity(this, val) - this.BlueSaturationIntensity = val; - end - function ret = get.BlueSaturationIntensity(this) - ret = this.BlueSaturationIntensity; - end - - function set.OrangePower(this,val) - this.OrangePower = val; - end - function ret = get.OrangePower(this) - ret = this.OrangePower; - end - function set.OrangeDetuning(this, val) - this.OrangeDetuning = val; - end - function ret = get.OrangeDetuning(this) - ret = this.OrangeDetuning; - end - function set.OrangeBeamRadius(this, val) - this.OrangeBeamRadius = val; - end - function ret = get.OrangeBeamRadius(this) - ret = this.OrangeBeamRadius; - end - function set.OrangeBeamWaist(this, val) - this.OrangeBeamWaist = val; - end - function ret = get.OrangeBeamWaist(this) - ret = this.OrangeBeamWaist; - end - function set.OrangeWaveNumber(this, val) - this.OrangeWaveNumber = val; - end - function ret = get.OrangeWaveNumber(this) - ret = this.OrangeWaveNumber; - end - function set.OrangeSaturationIntensity(this, val) - this.OrangeSaturationIntensity = val; - end - function ret = get.OrangeSaturationIntensity(this) - ret = this.OrangeSaturationIntensity; - end - - function set.CoolingBeamPower(this,val) - this.CoolingBeamPower = val; - end - function ret = get.CoolingBeamPower(this) - ret = this.CoolingBeamPower; - end - function set.CoolingBeamDetuning(this, val) - this.CoolingBeamDetuning = val; - end - function ret = get.CoolingBeamDetuning(this) - ret = this.CoolingBeamDetuning; - end - function set.CoolingBeamRadius(this, val) - this.CoolingBeamRadius = val; - end - function ret = get.CoolingBeamRadius(this) - ret = this.CoolingBeamRadius; - end - function set.CoolingBeamWaist(this, val) - this.CoolingBeamWaist = val; - end - function ret = get.CoolingBeamWaist(this) - ret = this.CoolingBeamWaist; - end - function set.CoolingBeamWaveNumber(this, val) - this.CoolingBeamWaveNumber = val; - end - function ret = get.CoolingBeamWaveNumber(this) - ret = this.CoolingBeamWaveNumber; - end - function set.CoolingBeamLinewidth(this, val) - this.CoolingBeamLinewidth = val; - end - function ret = get.CoolingBeamLinewidth(this) - ret = this.CoolingBeamLinewidth; - end - function set.CoolingBeamSaturationIntensity(this, val) - this.CoolingBeamSaturationIntensity = val; - end - function ret = get.CoolingBeamSaturationIntensity(this) - ret = this.CoolingBeamSaturationIntensity; - end - - function set.SidebandPower(this,val) - this.SidebandPower = val; - end - function ret = get.SidebandPower(this) - ret = this.SidebandPower; - end - function set.SidebandDetuning(this, val) - this.SidebandDetuning = val; - end - function ret = get.SidebandDetuning(this) - ret = this.SidebandDetuning; - end - function set.SidebandBeamRadius(this, val) - this.SidebandBeamRadius = val; - end - function ret = get.SidebandBeamRadius(this) - ret = this.SidebandBeamRadius; - end - function set.SidebandBeamWaist(this, val) - this.SidebandBeamWaist = val; - end - function ret = get.SidebandBeamWaist(this) - ret = this.SidebandBeamWaist; - end - function set.SidebandBeamSaturationIntensity(this, val) - this.SidebandBeamSaturationIntensity = val; - end - function ret = get.SidebandBeamSaturationIntensity(this) - ret = this.SidebandBeamSaturationIntensity; - end - - function set.PushBeamPower(this,val) - this.PushBeamPower = val; - end - function ret = get.PushBeamPower(this) - ret = this.PushBeamPower; - end - function set.PushBeamDetuning(this, val) - this.PushBeamDetuning = val; - end - function ret = get.PushBeamDetuning(this) - ret = this.PushBeamDetuning; - end - function set.PushBeamRadius(this, val) - this.PushBeamRadius = val; - end - function ret = get.PushBeamRadius(this) - ret = this.PushBeamRadius; - end - function set.PushBeamWaist(this, val) - this.PushBeamWaist = val; - end - function ret = get.PushBeamWaist(this) - ret = this.PushBeamWaist; - end - function set.PushBeamWaveNumber(this, val) - this.PushBeamWaveNumber= val; - end - function ret = get.PushBeamWaveNumber(this) - ret = this.PushBeamWaveNumber; - end - function set.PushBeamLinewidth(this, val) - this.PushBeamLinewidth = val; - end - function ret = get.PushBeamLinewidth(this) - ret = this.PushBeamLinewidth; - end - function set.PushBeamDistance(this, val) - this.PushBeamDistance = val; - end - function ret = get.PushBeamDistance(this) - ret = this.PushBeamDistance; - end - function set.DistanceBetweenPushBeamAnd3DMOTCenter(this, val) - this.DistanceBetweenPushBeamAnd3DMOTCenter = val; - end - function ret = get.DistanceBetweenPushBeamAnd3DMOTCenter(this) - ret = this.DistanceBetweenPushBeamAnd3DMOTCenter; - end - function set.PushBeamSaturationIntensity(this, val) - this.PushBeamSaturationIntensity = val; - end - function ret = get.PushBeamSaturationIntensity(this) - ret = this.PushBeamSaturationIntensity; - end - - function set.ZeemanSlowerBeamPower(this,val) - this.ZeemanSlowerBeamPower = val; - end - function ret = get.ZeemanSlowerBeamPower(this) - ret = this.ZeemanSlowerBeamPower; - end - function set.ZeemanSlowerBeamDetuning(this, val) - this.ZeemanSlowerBeamDetuning = val; - end - function ret = get.ZeemanSlowerBeamDetuning(this) - ret = this.ZeemanSlowerBeamDetuning; - end - function set.ZeemanSlowerBeamRadius(this, val) - this.ZeemanSlowerBeamRadius = val; - end - function ret = get.ZeemanSlowerBeamRadius(this) - ret = this.ZeemanSlowerBeamRadius; - end - function set.ZeemanSlowerBeamWaist(this, val) - this.ZeemanSlowerBeamWaist = val; - end - function ret = get.ZeemanSlowerBeamWaist(this) - ret = this.ZeemanSlowerBeamWaist; - end - function set.ZeemanSlowerBeamSaturationIntensity(this, val) - this.ZeemanSlowerBeamSaturationIntensity = val; - end - function ret = get.ZeemanSlowerBeamSaturationIntensity(this) - ret = this.ZeemanSlowerBeamSaturationIntensity; - end - - function set.TotalPower(this,val) - this.TotalPower = val; - end - function ret = get.TotalPower(this) - ret = this.TotalPower; - end - function set.LandegFactor(this,val) - this.LandegFactor = val; - end - function ret = get.LandegFactor(this) - ret = this.LandegFactor; - end - function set.MagneticSubLevel(this,val) - this.MagneticSubLevel = val; - end - function ret = get.MagneticSubLevel(this) - ret = this.MagneticSubLevel; - end - - function set.CaptureVelocity(this,val) - this.CaptureVelocity = val; - end - function ret = get.CaptureVelocity(this) - ret = this.CaptureVelocity; - end - function set.VelocityCutoff(this,val) - this.VelocityCutoff = val; - end - function ret = get.VelocityCutoff(this) - ret = this.VelocityCutoff; - end - function set.ClausingFactor(this,val) - this.ClausingFactor = val; - end - function ret = get.ClausingFactor(this) - ret = this.ClausingFactor; - end - function set.ReducedClausingFactor(this,val) - this.ReducedClausingFactor = val; - end - function ret = get.ReducedClausingFactor(this) - ret = this.ReducedClausingFactor; - end - function set.ReducedFlux(this,val) - this.ReducedFlux = val; - end - function ret = get.ReducedFlux(this) - ret = this.ReducedFlux; - end - function set.TimeSpentInInteractionRegion(this,val) - this.TimeSpentInInteractionRegion = val; - end - function ret = get.TimeSpentInInteractionRegion(this) - ret = this.TimeSpentInInteractionRegion; - end - - function set.DebugMode(this, val) - this.DebugMode = val; - end - function ret = get.DebugMode(this) - ret = this.DebugMode; - end - function set.DoSave(this, val) - this.DoSave = val; - end - 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 - function ret = get.CoolingBeamSaturationParameter(this) - ret = 0.1 * (4 * this.CoolingBeamPower) / (pi*this.CoolingBeamWaist^2 * this.CoolingBeamSaturationIntensity); % two beams are reflected - end - - function ret = get.SidebandSaturationParameter(this) - ret = 0.1 * (4 * this.SidebandPower) / (pi*this.SidebandBeamWaist^2 * this.SidebandBeamSaturationIntensity); - end - - function ret = get.PushBeamSaturationParameter(this) - ret = 0.1 * this.PushBeamPower/(pi * this.PushBeamWaist^2 * this.PushBeamSaturationIntensity); - end - - function ret = get.ZeemanSlowerBeamSaturationParameter(this) - ret = 0.1 * this.ZeemanSlowerBeamPower / (pi * this.ZeemanSlowerBeamWaist^2 * this.ZeemanSlowerBeamSaturationIntensity); - end - - function ret = get.OvenTemperatureinKelvin(this) - ret = this.OvenTemperature + Helper.PhysicsConstants.ZeroKelvin; - end - - function ret = get.AverageVelocity(this) - %See Background collision probability section in Barbiero - ret = sqrt((8 * pi *Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass)); - end - - function ret = get.AtomicBeamDensity(this) - %See Background collision probability section in Barbiero - ret = this.calculateFreeMolecularRegimeFlux() / (this.AverageVelocity * pi * (this.Beta*this.NozzleLength/2)^2); - end - - function ret = get.MeanFreePath(this) - % Cross section = pi ( 2 * Van-der-waals radius of Dy)^2; - % Van-der-waals radius of Dy = 281e-12 - %See Expected atomic flux section and Background collision probability section in Barbiero - ret = 1/(sqrt(2) * ( pi * (2*281e-12)^2) * this.AtomicBeamDensity); - end - - function ret = get.CollisionTime(this) - ret = 3 * this.MeanFreePath/this.AverageVelocity; %See Background collision probability section in Barbiero - end - - end % - getters for dependent properties - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - methods(Access = protected) - function cp = copyElement(this) - % Shallow copy object - cp = copyElement@matlab.mixin.Copyable(this); - - % Forces the setter to redefine the function handles to the new copied object - % cp.potentialType = this.potentialType; - - pl = properties(this); - for k = 1:length(pl) - sc = superclasses(this.(pl{k})); - if any(contains(sc,{'matlab.mixin.Copyable'})) - cp.(pl{k}) = this.(pl{k}).copy(); - end - end - end - end - - methods (Static) - - % Creates an Instance of Class, ensures singleton behaviour (that there - % can only be one Instance of this class - function singleObj = getInstance(varargin) - % Creates an Instance of Class, ensures singleton behaviour - persistent localObj; - if isempty(localObj) || ~isvalid(localObj) - localObj = MOTSimulator(varargin{:}); - end - singleObj = localObj; - end - end - -end diff --git a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m deleted file mode 100644 index c3ce154..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m +++ /dev/null @@ -1,27 +0,0 @@ -function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector) - % is the distance between the chamber center and the cross point of push beam and z-axis (along the gravity) - WaveVectorEndPoint = [0, 1, this.DistanceBetweenPushBeamAnd3DMOTCenter/this.PushBeamDistance]; - WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint); - Origin=[0,0,0]; - - SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist); - - DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * this.PushBeamWaveNumber; - - Detuning = this.PushBeamDetuning - DopplerShift; - - s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2)); - a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5); - a_push = a_sat .* (s_push/(1+s_push)); - - if this.SpontaneousEmission - a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, this.PushBeamLinewidth, this.PushBeamWaveNumber); - else - a_scatter = [0,0,0]; - end - - a_total = a_push + a_scatter; - - ret = a_total(1:3); -end - diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m b/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m deleted file mode 100644 index f235207..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m +++ /dev/null @@ -1,28 +0,0 @@ -function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector) - - switch this.SimulationMode - case "2D" - VelocityUnitVector = VelocityVector./norm(VelocityVector); - UpperLimit = 500; - LowerLimit = 0; - - for Index = 1:500 - InitialVelocity = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; - ParticleDynamicalQuantities = this.solver(PositionVector, InitialVelocity); - FinalPositionVector = ParticleDynamicalQuantities(end, 1:3); - if rssq(FinalPositionVector) <= this.OvenDistance - LowerLimit = 0.5 * (UpperLimit + LowerLimit); - else - UpperLimit = 0.5 * (UpperLimit + LowerLimit); - end - - if UpperLimit - LowerLimit < 1 - ret = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; - break; - end - end - case "3D" - % Development In progress - end -end - diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m deleted file mode 100644 index 89ba894..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m +++ /dev/null @@ -1,39 +0,0 @@ -function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this) - switch this.SimulationMode - case "2D" - n = this.NumberOfAtoms; - DynamicalQuantities = this.ParticleDynamicalQuantities; - NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep); - NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps); - TimeCounts = zeros(1, n); - CollisionEvents = zeros(1, n); - - % Include the stochastic process of background collisions - for AtomIndex = 1:n - TimeCounts(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(DynamicalQuantities(AtomIndex,:,1:3))); - end - this.TimeSpentInInteractionRegion = mean(TimeCounts); - for AtomIndex = 1:n - CollisionEvents(AtomIndex) = this.computeCollisionProbability(); - end - - % Count the number of loaded atoms subject to conditions - for TimeIndex = 1:NumberOfTimeSteps - if TimeIndex ~= 1 - NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1); - end - for AtomIndex = 1:n - Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))'; - Velocity = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 4:6))'; - if this.exitCondition(Position, Velocity, CollisionEvents(AtomIndex)) - NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1; - end - end - end - - [LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(NumberOfLoadedAtoms); - - case "3D" - % Development In progress - end -end diff --git a/MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m b/MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m deleted file mode 100644 index 9bf974d..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m +++ /dev/null @@ -1,67 +0,0 @@ -function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = 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); - ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2); - - 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, ConfidenceInterval] = this.runSimulation(); - LoadingRateArray(i) = LoadingRate; - StandardErrorArray(i) = StandardError; - ConfidenceIntervalArray(i,1) = ConfidenceInterval(1); - ConfidenceIntervalArray(i,2) = ConfidenceInterval(2); - end - - if this.DoSave - LoadingRate = struct; - LoadingRate.Values = LoadingRateArray; - LoadingRate.Errors = StandardErrorArray; - LoadingRate.CI = ConfidenceIntervalArray; - 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 deleted file mode 100644 index a930592..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m +++ /dev/null @@ -1,82 +0,0 @@ -function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = 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); - - addParameter(p, 'ChangeInitialConditions', false, @islogical); - addParameter(p, 'ParameterNameArray', {}, @iscell); - addParameter(p, 'ParameterValueArray', {}, @iscell); - - 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; - ChangeInitialConditions = p.Results.ChangeInitialConditions; - ParameterNameArray = p.Results.ParameterNameArray; - ParameterValueArray = p.Results.ParameterValueArray; - - 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('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 - 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, StandardError, ConfidenceInterval] = this.runSimulation(); - LoadingRateArray(i, j) = LoadingRate; - StandardErrorArray(i, j) = StandardError; - ConfidenceIntervalArray(i, j, 1) = ConfidenceInterval(1); - ConfidenceIntervalArray(i, j, 2) = ConfidenceInterval(2); - end - end - - if this.DoSave - LoadingRate = struct; - LoadingRate.Values = LoadingRateArray; - LoadingRate.Errors = StandardErrorArray; - LoadingRate.CI = ConfidenceIntervalArray; - 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/initialPositionSampling.m b/MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m deleted file mode 100644 index 197912f..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m +++ /dev/null @@ -1,12 +0,0 @@ -function ret = initialPositionSampling(this) - switch this.SimulationMode - case "2D" - n = this.NumberOfAtoms; - phi = 2 * pi * rand(n,1); - rho = this.Beta * 0.5 * this.NozzleLength * sqrt(rand(n,1)); - ret = [-this.OvenDistance * ones(n,1), rho.*cos(phi), rho.*sin(phi)]; - case "3D" - % Development In progress - end -end - diff --git a/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m b/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m deleted file mode 100644 index 7235e1f..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m +++ /dev/null @@ -1,67 +0,0 @@ -function ret = initialVelocitySampling(this) - switch this.SimulationMode - case "2D" - n = this.NumberOfAtoms; - - % Calculate Calculate Capture velocity --> Introduce velocity cutoff - - this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); - 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. - - [ReducedClausingFactor, NormalizationConstantForAngularDistribution] = this.calculateReducedClausingFactor(); - this.ReducedClausingFactor = ReducedClausingFactor; - - VelocityDistribution = @(velocity) sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin))^3 ... - * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ... - * this.OvenTemperatureinKelvin))); - - c = integral(VelocityDistribution, 0, this.VelocityCutoff); - - this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux(); - - ret = zeros(n,3); - 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 > this.VelocityCutoff - MaximumVelocityAllowed = this.VelocityCutoff; - else - MaximumVelocityAllowed = MostProbableVelocity; - end - NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed); - - parfor i = 1:n - % Rejection Sampling of speed - y = rand(1); - x = this.VelocityCutoff * rand(1); - while y > ((NormalizationConstantForVelocityDistribution)^-1 * this.velocityDistributionFunction(x)) %As long as this loop condition is satisfied, reject the corresponding x value - y = 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 = rand(1); - z = this.NozzleExitDivergence * rand(1); - - while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value - w = rand(1); - z = this.NozzleExitDivergence * rand(1); - end - SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample - - % Sampling of azimuthal angle - SampledAzimuthalAngle(i)= 2 * pi * rand(1); - - ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ... - SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))]; - end - case "3D" - % Development In progress - end -end - diff --git a/MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m b/MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m deleted file mode 100644 index 802c27b..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m +++ /dev/null @@ -1,14 +0,0 @@ -function ret = magneticFieldForMOT(this, r) - switch this.SimulationMode - case '2D' - ret = zeros(1,4); - alpha = this.MagneticGradient; - ret(1) = r(3)*alpha; - ret(2) = 0; - ret(3) = r(1)*alpha; - ret(4) = sqrt(ret(1)^2+ret(2)^2+ret(3)^2); - case '3D' - % Development in progress - - end -end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m deleted file mode 100644 index 71a6544..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m +++ /dev/null @@ -1,35 +0,0 @@ -function reinitializeSimulator(this) - %% PHYSICAL CONSTANTS - pc = Helper.PhysicsConstants; - %% SIMULATION PARAMETERS - this.NozzleLength = 60e-3; - this.NozzleRadius = 2.60e-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); - % 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 = 0.32; % Distance between the 2-D MOT the 3-D MOT - this.BlueWaveNumber = 2*pi/pc.BlueWavelength; - this.BlueSaturationIntensity = 0.1 * (2 * pi^2 / 3) * ((pc.PlanckConstantReduced * pc.SpeedOfLight * pc.BlueLinewidth) / (pc.BlueWavelength)^3); - this.OrangeWaveNumber = 2*pi/pc.OrangeWavelength; - this.OrangeSaturationIntensity = 0.1 * (2 * pi^2 / 3) * ((pc.PlanckConstantReduced * pc.SpeedOfLight * pc.OrangeLinewidth) / (pc.OrangeWavelength)^3); - this.BlueBeamRadius = min(0.035/2,sqrt(2)/2*this.OvenDistance); % Diameter of CF40 flange = 0.035 - Theta_Nozzle = atan((this.NozzleRadius+this.BlueBeamRadius*sqrt(2))/this.OvenDistance); % The angle of capture region towards the oven nozzle - Theta_Aperture = 15/360*2*pi; % The limitation angle of the second aperture in the oven - this.NozzleExitDivergence = min(Theta_Nozzle,Theta_Aperture); - this.MOTExitDivergence = 16e-3; % The limitation angle between 2D-MOT and 3D-MOT - this.TotalPower = 0.6; - this.OrangeBeamRadius = 1.2e-3; - this.PushBeamRadius = 1.2e-3; - this.PushBeamDistance = 0.32; - this.PushBeamLinewidth = Helper.PhysicsConstants.OrangeLinewidth; - this.PushBeamWaveNumber = this.OrangeWaveNumber; - this.PushBeamSaturationIntensity = this.OrangeSaturationIntensity; - this.ZeemanSlowerBeamRadius = 0.035; - this.ZeemanSlowerBeamSaturationIntensity = this.BlueSaturationIntensity; - this.DistanceBetweenPushBeamAnd3DMOTCenter = 0; -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 deleted file mode 100644 index b18dd3f..0000000 --- a/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m +++ /dev/null @@ -1,137 +0,0 @@ -function setInitialConditions(this, varargin) - - p = inputParser; - p.KeepUnmatched = true; - addParameter(p, 'NumberOfAtoms', 5000,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'BluePower', this.TotalPower,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'BlueDetuning', -1.39*Helper.PhysicsConstants.BlueLinewidth,... - @(x) assert(isnumeric(x) && isscalar(x))); - addParameter(p, 'BlueBeamWaist', 16.6667e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'SidebandPower', 0.5*this.TotalPower,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'SidebandDetuning', -4.5*Helper.PhysicsConstants.BlueLinewidth,... - @(x) assert(isnumeric(x) && isscalar(x))); - addParameter(p, 'SidebandBeamWaist', 15e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'PushBeamPower', 100e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'PushBeamDetuning', -5*Helper.PhysicsConstants.OrangeLinewidth,... - @(x) assert(isnumeric(x) && isscalar(x))); - addParameter(p, 'PushBeamWaist', 0.81e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'OrangePower', 70e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'OrangeDetuning', -1*Helper.PhysicsConstants.OrangeLinewidth,... - @(x) assert(isnumeric(x) && isscalar(x))); - addParameter(p, 'OrangeBeamWaist', 12e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'ZeemanSlowerBeamPower', 200e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'ZeemanSlowerBeamDetuning', -7*Helper.PhysicsConstants.BlueLinewidth,... - @(x) assert(isnumeric(x) && isscalar(x))); - addParameter(p, 'ZeemanSlowerBeamWaist', 7e-3,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - 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; - this.BlueDetuning = p.Results.BlueDetuning; - this.BlueBeamWaist = p.Results.BlueBeamWaist; - this.SidebandPower = p.Results.SidebandPower; - this.SidebandDetuning = p.Results.SidebandDetuning; - this.SidebandBeamWaist = p.Results.SidebandBeamWaist; - this.PushBeamPower = p.Results.PushBeamPower; - this.PushBeamDetuning = p.Results.PushBeamDetuning; - this.PushBeamWaist = p.Results.PushBeamWaist; - this.OrangePower = p.Results.OrangePower; - this.OrangeDetuning = p.Results.OrangeDetuning; - this.OrangeBeamWaist = p.Results.OrangeBeamWaist; - this.ZeemanSlowerBeamPower = p.Results.ZeemanSlowerBeamPower; - this.ZeemanSlowerBeamDetuning = p.Results.ZeemanSlowerBeamDetuning; - this.ZeemanSlowerBeamWaist = p.Results.ZeemanSlowerBeamWaist; - this.MagneticGradient = p.Results.MagneticGradient; - %% Set general parameters according to simulation mode - switch this.SimulationMode - case "2D" - this.CoolingBeamPower = this.BluePower; - this.CoolingBeamWaist = this.BlueBeamWaist; - this.CoolingBeamLinewidth = Helper.PhysicsConstants.BlueLinewidth; - this.CoolingBeamWaveNumber = this.BlueWaveNumber; - this.CoolingBeamDetuning = this.BlueDetuning; - this.CoolingBeamRadius = this.BlueBeamRadius; - this.CoolingBeamWaist = this.BlueBeamWaist; - this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity; - this.SidebandBeamRadius = this.BlueBeamRadius; - this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity; - this.LandegFactor = Helper.PhysicsConstants.BlueLandegFactor; - this.MagneticSubLevel = 1; - case "3D" - % Development In progress - end - - - %% - store in struct - 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/test_MOTSimulator.m b/MOT Capture Process Simulation/test_MOTSimulator.m deleted file mode 100644 index 31a690a..0000000 --- a/MOT Capture Process Simulation/test_MOTSimulator.m +++ /dev/null @@ -1,188 +0,0 @@ -%% - Create solver object with specified options -clc -%% - -OptionsStruct = struct; -OptionsStruct.SimulationMode = '2D'; -OptionsStruct.TimeStep = 50e-06; % in s -OptionsStruct.SimulationTime = 4e-03; % in s -OptionsStruct.SpontaneousEmission = true; -OptionsStruct.Sideband = false; -OptionsStruct.PushBeam = false; -OptionsStruct.Gravity = false; -OptionsStruct.BackgroundCollision = true; -OptionsStruct.SaveData = false; -OptionsStruct.SaveDirectory = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation'; -options = Helper.convertstruct2cell(OptionsStruct); - -Simulator = MOTSimulator(options{:}); - -clear OptionsStruct -%% - Set Initial Conditions: Run with default values -Simulator.setInitialConditions(); - -%% - Set Initial Conditions: Set manually -OptionsStruct = struct; -OptionsStruct.NumberOfAtoms = 5000; -OptionsStruct.BluePower = 0.2; % in W -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 -OptionsStruct.PushBeamPower = 0.010; % in W -OptionsStruct.PushBeamDetuning = 0; % in Hz -OptionsStruct.PushBeamWaist = 0.005; % in m - -options = Helper.convertstruct2cell(OptionsStruct); -Simulator.setInitialConditions(options{:}); -clear OptionsStruct -%% - Run Simulation -[LoadingRate, ~] = Simulator.runSimulation(); - -%% - Plot initial distribution -Simulator.setInitialConditions(); -% - sampling the position distribution -InitialPositions = Simulator.initialPositionSampling(); -% - sampling the velocity distribution -InitialVelocities = Simulator.initialVelocitySampling(); -NumberOfBins = 100; -Plotting.plotPositionAndVelocitySampling(NumberOfBins, InitialPositions, InitialVelocities); - -%% - Plot distributions of magnitude and direction of initial velocities -NumberOfBins = 50; -Plotting.plotInitialVeloctiySamplingVsAngle(Simulator, NumberOfBins) - -%% - Plot Magnetic Field -XAxisRange = [-5 5]; -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.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius) - -%% - Plot the Free Molecular Flux for different temperatures -Temperature = [950, 1000, 1050]; % 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(); -Plotting.plotCaptureVelocityVsAngle(Simulator); - -%% - Plot Phase Space - -Simulator.NumberOfAtoms = 100; -MinimumVelocity = 0; -MaximumVelocity = 150; -NumberOfBins = 200; %Along each axis -IncidentAtomDirection = 0*2*pi/360; -IncidentAtomPosition = 0; -Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition) - -%% - Plot trajectories along the 3 directions - -Simulator.NumberOfAtoms = 100; -MaximumVelocity = 150; -IncidentAtomDirection = 0*2*pi/360; -IncidentAtomPosition = 0; -Plotting.plotDynamicalQuantities(Simulator, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition); -%% - Scan parameters - -% ONE-PARAMETER SCAN - -NumberOfPointsForParam = 5; %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; - -OptionsStruct = struct; -OptionsStruct.ChangeInitialConditions = true; -OptionsStruct.ParameterNameArray = {'NumberOfAtoms'}; -OptionsStruct.ParameterValueArray = {5000}; -options = Helper.convertstruct2cell(OptionsStruct); - - -tStart = tic; -[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:}); -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)'; -OptionsStruct.RescalingFactorForYQuantity = 1e-11; -OptionsStruct.ErrorsForYQuantity = true; -OptionsStruct.ErrorsArray = StandardErrorArray; -OptionsStruct.CIForYQuantity = true; -OptionsStruct.CIArray = ConfidenceIntervalArray; -OptionsStruct.RemoveOutliers = true; -OptionsStruct.YLabelString = 'Loading rate (x 10^{11} 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 = 10; %iterations of the simulation -NumberOfPointsForSecondParam = 10; - -% Scan Sideband Detuning and Power Ratio -DetuningArray = linspace(-0.5,-10, NumberOfPointsForParam) * 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; -OptionsStruct.ChangeInitialConditions = true; -OptionsStruct.ParameterNameArray = {'NumberOfAtoms'}; -OptionsStruct.ParameterValueArray = {5000}; -options = Helper.convertstruct2cell(OptionsStruct); - -tStart = tic; -[LoadingRateArray, StandardErrorArray, ConfidenceInterval] = 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; -SecondParameterArray = SidebandPowerArray; -QuantityOfInterestArray = LoadingRateArray; - -OptionsStruct = struct; -OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1; -OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)'; -OptionsStruct.RescalingFactorForSecondParameter = 1000; -OptionsStruct.YLabelString = 'Sideband Power (mW)'; -OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-10; -OptionsStruct.ZLabelString = 'Loading rate (x 10^{10} atoms/s)'; -OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100); - -options = Helper.convertstruct2cell(OptionsStruct); - -Plotting.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, QuantityOfInterestArray, options{:}) - -clear OptionsStruct \ No newline at end of file diff --git a/test_MOTCaptureProcessSimulation.m b/test_MOTCaptureProcessSimulation.m new file mode 100644 index 0000000..95db29e --- /dev/null +++ b/test_MOTCaptureProcessSimulation.m @@ -0,0 +1,188 @@ +%% This script is testing the functionalities of the MOT Capture Process Simulation Classes +% +% Important: Run only sectionwise!! + +%% - Testing the MOTCaptureProcess-Class +% - Create MOTCaptureProcess object with specified options +% - Automatically creates Beams objects +OptionsStruct = struct; +OptionsStruct.NumberOfAtoms = 10000; +OptionsStruct.TimeStep = 50e-06; % in s +OptionsStruct.SimulationTime = 4e-03; % in s +OptionsStruct.SpontaneousEmission = true; +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'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +Oven = Simulator.Oven(options{:}); +MOT2D = Simulator.TwoDimensionalMOT(options{:}); +Beams = MOT2D.Beams; + +%% - Run Simulation +poolobj = gcp('nocreate'); % Check if pool is open +if isempty(poolobj) + parpool; +end +[LoadingRate, ~] = MOT2D.runSimulation(Oven); +%% - Plot initial distribution +% - sampling the position distribution +InitialPositions = Oven.initialPositionSampling(); +% - sampling the velocity distribution +InitialVelocities = Oven.initialVelocitySampling(MOT2D); +NumberOfBins = 100; +Plotter.plotPositionAndVelocitySampling(NumberOfBins, InitialPositions, InitialVelocities); + +%% - Plot distributions of magnitude and direction of initial velocities +NumberOfBins = 50; +Plotter.plotInitialVeloctiySamplingVsAngle(Oven, MOT2D, NumberOfBins) + +%% - Plot Magnetic Field +XAxisRange = [-5 5]; +YAxisRange = [-5 5]; +ZAxisRange = [-5 5]; +Plotter.visualizeMagneticField(MOT2D, XAxisRange, YAxisRange, ZAxisRange) + +%% - Plot MFP & VP for different temperatures +TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius +Plotter.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius) + +%% - Plot the Free Molecular Flux for different temperatures +Temperature = [950, 1000, 1050]; % Temperature +Plotter.plotFreeMolecularFluxVsTemp(Oven,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 +Plotter.plotAngularDistributionForDifferentBeta(Oven, Beta) + +%% - Plot Capture Velocity +Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot! + +%% - Plot Phase Space with Acceleration Field +MOT2D.Sideband = false; +MOT2D.NumberOfAtoms = 50; +MinimumVelocity = 0; +MaximumVelocity = 150; +NumberOfBins = 200; %Along each axis +IncidentAtomDirection = 0*2*pi/360; +IncidentAtomPosition = 0; +Plotter.plotPhaseSpaceWithAccelerationField(Oven, MOT2D, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition) + +%% - Plot Trajectories along the 3 directions +MOT2D.NumberOfAtoms = 100; +MaximumVelocity = 150; +IncidentAtomDirection = 0*2*pi/360; +IncidentAtomPosition = 0; + +%% - Positions +Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, 'PlotPositions', true); + +%% - Velocities +Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, 'PlotVelocities', true); + +%% - Scan parameters: One-Parameter Scan + +MOT2D.NumberOfAtoms = 5000; +MOT2D.TotalPower = 0.4; +CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)}; + +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); + +tStart = tic; +for i=1:NumberOfPointsForFirstParam + CoolingBeam.Power = PowerArray(i); + [LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i, :)] = MOT2D.runSimulation(Oven); +end +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)'; +OptionsStruct.RescalingFactorForYQuantity = 1e-10; +OptionsStruct.ErrorsForYQuantity = true; +OptionsStruct.ErrorsArray = StandardErrorArray; +OptionsStruct.CIForYQuantity = true; +OptionsStruct.CIArray = ConfidenceIntervalArray; +OptionsStruct.RemoveOutliers = true; +OptionsStruct.YLabelString = 'Loading rate (x 10^{10} atoms/s)'; +OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100); + +options = Helper.convertstruct2cell(OptionsStruct); + +Plotter.plotResultForOneParameterScan(ParameterArray, QuantityOfInterestArray, options{:}) + +clear OptionsStruct + +%% - Scan parameters: Two-Parameter Scan + +MOT2D.NumberOfAtoms = 5000; +MOT2D.TotalPower = 0.6; +MOT2D.Sideband = true; +SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)}; + +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; + +LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); +StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); +ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2); + +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 +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.RescalingFactorForSecondParameter = 1000; +OptionsStruct.YLabelString = 'Sideband Power (mW)'; +OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-10; +OptionsStruct.ZLabelString = 'Loading rate (x 10^{10} atoms/s)'; +OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100); + +options = Helper.convertstruct2cell(OptionsStruct); + +Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, QuantityOfInterestArray, options{:}) + +clear OptionsStruct \ No newline at end of file