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

Compare commits

base: karthik:c2e6420243e31a3385f4b9a311aa7ca60f0eae9c
Branches Tags
karthik:master
..
compare: karthik:d0b4d0b9b84f140b55cb75815a6aa6c983b51a4f
Branches Tags
karthik:master

No commits in common. "c2e6420243e31a3385f4b9a311aa7ca60f0eae9c" and "d0b4d0b9b84f140b55cb75815a6aa6c983b51a4f" have entirely different histories.
c2e6420243 ... d0b4d0b9b8

35 changed files with 403 additions and 601 deletions
Show Stats Expand all files Collapse all files

10
MOT Capture Process Simulation/+Helper/calculateDistanceFromPointToLine.m
View File

@ -1,10 +0,0 @@
function ret = calculateDistanceFromPointToLine(p0 , p1, p2)
p01 = p0 - p1;
p12 = p2 - p1;
CrossProduct = [p01(2)*p12(3) - p01(3)*p12(2), p01(3)*p12(1) - p01(1)*p12(3), p01(1)*p12(2) - p01(2)*p12(1)];
ret = rssq(CrossProduct) / rssq(p12);
%Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
%Area = One side of parallelogram X Base
%ret = norm(cross(one side, base))./ norm(base);
end

18
MOT Capture Process Simulation/+Helper/findAllZeroCrossings.m
View File

@ -1,18 +0,0 @@
function ret = findAllZeroCrossings(x,y)
% Finds all Zero-crossing of the function y = f(x)
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
zxidx = zci(y);
if ~isempty(zxidx)
for k1 = 1:numel(zxidx)
idxrng = max([1 zxidx(k1)-1]):min([zxidx(k1)+1 numel(y)]);
xrng = x(idxrng);
yrng = y(idxrng);
[yrng2, ~, jyrng] = unique(yrng); %yrng is a new array containing the unique values of yrng. jyrng contains the indices in yrng that correspond to the original vector. yrng = yrng2(jyrng)
xrng2 = accumarray(jyrng, xrng, [], @mean); %This function creates a new array "xrng2" by applying the function "@mean" to all elements in "xrng" that have identical indices in "jyrng". Any elements with identical X values will have identical indices in jyrng. Thus, this function creates a new array by averaging values with identical X values in the original array.
ret(k1) = interp1( yrng2(:), xrng2(:), 0, 'linear', 'extrap' );
end
else
warning('No zero crossings found!')
ret = nan;
end
end

2
MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m
View File

@ -64,7 +64,7 @@ function plotAngularDistributionForDifferentBeta(obj, Beta)
leg = legend; leg = legend;
hXLabel = xlabel('\theta (rad)'); hXLabel = xlabel('\theta (rad)');
hYLabel = ylabel('J(\theta)'); hYLabel = ylabel('J(\theta)');
hTitle = sgtitle('Angular Distribution (Transition Flow)'); hTitle = sgtitle('Angular Distribution (Free Molecular Flow)');
set([hXLabel, hYLabel, leg] , ... set([hXLabel, hYLabel, leg] , ...
'FontSize' , 14 ); 'FontSize' , 14 );

10
MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m
View File

@ -1,6 +1,6 @@
function plotCaptureVelocityVsAngle(obj) function plotCaptureVelocityVsAngle(obj)
theta = linspace(0, obj.NozzleExitDivergence, 1000); theta = linspace(0, obj.MOTExitDivergence, 40);
CaptureVelocity = zeros(length(theta),3); CaptureVelocity = zeros(length(theta),3);
for i=1:length(theta) for i=1:length(theta)
@ -20,16 +20,16 @@ function plotCaptureVelocityVsAngle(obj)
screensize = get(0,'ScreenSize'); screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600]; f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
plot(theta .* 1e+03, sqrt(CaptureVelocity(:,1).^2+CaptureVelocity(:,2).^2+CaptureVelocity(:,3).^2), 'Linewidth', 1.5) plot(theta, sqrt(CaptureVelocity(:,1).^2+CaptureVelocity(:,2).^2+CaptureVelocity(:,3).^2), 'Linewidth', 1.5)
hXLabel = xlabel('\theta (mrad)'); hXLabel = xlabel('\theta (rad)');
hYLabel = ylabel('Velocity (m/s)'); hYLabel = ylabel('Velocity (m/s)');
hTitle = sgtitle('Capture velocity for different angles of divergence'); hTitle = sgtitle('Capture Velocity of an atomic beam from (0,0,0)');
set([hXLabel, hYLabel] , ... set([hXLabel, hYLabel] , ...
'FontSize' , 14 ); 'FontSize' , 14 );
set( hTitle , ... set( hTitle , ...
'FontSize' , 14 ); 'FontSize' , 18 );
grid on grid on
Helper.bringFiguresWithTagInForeground(); Helper.bringFiguresWithTagInForeground();

18
MOT Capture Process Simulation/+Plotting/plotConfidenceIntervalRegion.m
View File

@ -1,18 +0,0 @@
function plotConfidenceIntervalRegion(x, y1, y2)
% draws two lines on a plot and shades the area between those
% lines to indicate confidence interval.
hold on
X_interpolated = linspace(min(x), max(x), 100);
Y1_interpolated = interp1(x,y1,X_interpolated);
Y2_interpolated = interp1(x,y2,X_interpolated);
%Plot the line edges
%plot(X_interpolated, Y1_interpolated, 'LineWidth', 0.5, 'LineStyle', '--', 'Color', '#FE1A1A');
%plot(X_interpolated, Y2_interpolated, 'LineWidth', 0.5, 'LineStyle', '--', 'Color', '#FE1A1A');
fill([X_interpolated fliplr(X_interpolated)], [Y1_interpolated fliplr(Y2_interpolated)], [0 71 138] ./ 255, 'EdgeColor', 'none', 'FaceAlpha', 0.2);
hold off
end

13
MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m
View File

@ -17,22 +17,21 @@ function plotFreeMolecularFluxVsTemp(obj, Temperature)
for i=1:length(Temperature) for i=1:length(Temperature)
beta = linspace(0.01,0.5,200); beta = linspace(0.01,0.5,200);
L = 2*obj.NozzleRadius./beta;
obj.OvenTemperature = Temperature(i); obj.OvenTemperature = Temperature(i);
flux = zeros(1,length(beta)); flux = zeros(1,length(L));
for j=1:length(beta) for j=1:length(L)
obj.Beta = beta(j); obj.NozzleLength = L(j);
[ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor(); flux(j) = obj.calculateFreeMolecularRegimeFlux();
flux(j) = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux();
end end
plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5) plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5)
end end
set(gca,'yscale','log') set(gca,'yscale','log')
obj.reinitializeSimulator(); obj.reinitializeSimulator();
[ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor();
xline(obj.Beta, 'k--','Linewidth', 0.5); xline(obj.Beta, 'k--','Linewidth', 0.5);
fmf = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux(); fmf = obj.calculateFreeMolecularRegimeFlux();
yline(fmf, 'k--', 'Linewidth', 1.5); yline(fmf, 'k--', 'Linewidth', 1.5);
textstring = [sprintf('%1.e',fmf) ' atoms/s for ' '$$ \beta $$ = ' num2str(obj.Beta, '%.2f') sprintf(' @ %.2f K', obj.OvenTemperatureinKelvin)]; 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'); txt = text((obj.Beta + 0.05*obj.Beta), (max(fmf) + 0.2*fmf), textstring, 'Interpreter','latex');

74
MOT Capture Process Simulation/+Plotting/plotInitialVeloctiySamplingVsAngle.m
View File

@ -1,74 +0,0 @@
function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins)
n = obj.NumberOfAtoms;
obj.setInitialConditions();
VelocitySamples = initialVelocitySampling(obj);
Speeds = zeros(length(VelocitySamples),1);
Angles = zeros(length(VelocitySamples),1);
for i=1:length(VelocitySamples)
Speeds(i) = sqrt(VelocitySamples(i,1)^2+VelocitySamples(i,2)^2+VelocitySamples(i,3)^2);
Angles(i) = acos(VelocitySamples(i,1)/Speeds(i));
end
SpeedsArray = linspace(0,obj.VelocityCutoff,5000);
SpeedBinSize = obj.VelocityCutoff / NumberOfBins;
VelocityDistribution = @(velocity) sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * obj.OvenTemperatureinKelvin))^3 ...
* velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* obj.OvenTemperatureinKelvin)));
c = integral(VelocityDistribution, 0, obj.VelocityCutoff);
AnglesArray = linspace(0, obj.NozzleExitDivergence,1000);
AngleBinSize = obj.NozzleExitDivergence / NumberOfBins;
dtheta = AnglesArray(2)-AnglesArray(1);
AngularDistribution = zeros(1,length(AnglesArray));
d = 0;
for i = 1:length(AnglesArray)
AngularDistribution(i) = obj.angularDistributionFunction(AnglesArray(i));
d = d + (2 * pi * AngularDistribution(i) * sin(AnglesArray(i)) * dtheta);
end
AngularDistribution = AngularDistribution .* (2 * pi .* sin(AnglesArray));
f_h = Helper.getFigureByTag('Velocity Sampling');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Velocity Sampling';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/10 screensize(4)/4] 1570,600];
subplot(1,2,1)
h_1 = histogram(Speeds, NumberOfBins,'FaceAlpha',0.4, 'EdgeAlpha',0.4);
hold on
plot(SpeedsArray, n/c * SpeedBinSize .* VelocityDistribution(SpeedsArray),'--', 'Linewidth',1.5)
xline(obj.VelocityCutoff, 'k--', 'Linewidth', 1.5);
text((obj.VelocityCutoff - 0.2 * obj.VelocityCutoff), max(h_1.Values) + 0.01 * max(h_1.Values), 'Cut-Off ---------->');
hXLabel_1 = xlabel('Magnitudes (m/s)');
hYLabel_1 = ylabel('Counts');
hold off
subplot(1,2,2)
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 ----------->');
hXLabel_2 = xlabel(' Directions (mrad)');
hYLabel_2 = ylabel('Counts');
hold off
hTitle = sgtitle('Sampled Velocities');
set([hXLabel_1, hXLabel_2, hYLabel_1, hYLabel_2] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

6
MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m
View File

@ -50,20 +50,20 @@ function plotPhaseSpaceWithAccelerationField(obj, MaximumVelocity, NumberOfBins,
%------------------------------------------------------------------------- %-------------------------------------------------------------------------
Y = linspace(0,MaximumVelocity,N); Y = linspace(0,MaximumVelocity,N);
ParticleDynamicalQuantities = zeros(length(Y),int64(T/tau),6); Trajectories = zeros(length(Y),int64(T/tau),9);
for i=1:length(Y) for i=1:length(Y)
x =-L/2; x =-L/2;
vx = Y(i)*cos(Theta); vx = Y(i)*cos(Theta);
vz = Y(i)*sin(Theta); vz = Y(i)*sin(Theta);
r = [x,0,z]; r = [x,0,z];
v = [vx,0,vz]; v = [vx,0,vz];
ParticleDynamicalQuantities(i,:,:) = obj.solver(r, v); [Trajectories(i,:,:),~] = obj.solver(r, v);
end end
hold on hold on
for i=1:length(Y) for i=1:length(Y)
plot(ParticleDynamicalQuantities(i,:,1),ParticleDynamicalQuantities(i,:,4),'w','linewidth',1.3) plot(Trajectories(i,:,1),Trajectories(i,:,4),'w','linewidth',1.3)
end end
hold off hold off

8
MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m
View File

@ -41,14 +41,14 @@ function plotPositionAndVelocitySampling(obj, NumberOfBins)
legend('FontSize', 14) legend('FontSize', 14)
subplot(3,2,4) subplot(3,2,4)
histogram(initialVelocities(:, 2),NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component') histogram(initialVelocities(:, 2)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component')
xlabel('Velocities (m/s)','FontSize', 14) xlabel('Velocities (mm/s)','FontSize', 14)
ylabel('Counts','FontSize', 14) ylabel('Counts','FontSize', 14)
legend('FontSize', 14) legend('FontSize', 14)
subplot(3,2,6) subplot(3,2,6)
histogram(initialVelocities(:, 3),NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component') histogram(initialVelocities(:, 3)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component')
xlabel('Velocities (m/s)','FontSize', 14) xlabel('Velocities (mm/s)','FontSize', 14)
ylabel('Counts','FontSize', 14) ylabel('Counts','FontSize', 14)
legend('FontSize', 14) legend('FontSize', 14)

30
MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m
View File

@ -10,10 +10,7 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
addParameter(p, 'XLabelString', 'X parameter', @ischar) addParameter(p, 'XLabelString', 'X parameter', @ischar)
addParameter(p, 'ErrorsForYQuantity', false, @islogical) addParameter(p, 'ErrorsForYQuantity', false, @islogical)
addParameter(p, 'ErrorsArray', [], @isvector) addParameter(p, 'ErrorsArray', [], @isvector)
addParameter(p, 'CIForYQuantity', false, @islogical)
addParameter(p, 'CIArray', [], @ismatrix)
addParameter(p, 'RescalingFactorForYQuantity', 1, @isscalar) addParameter(p, 'RescalingFactorForYQuantity', 1, @isscalar)
addParameter(p, 'RemoveOutliers', false, @islogical)
addParameter(p, 'YLabelString', 'Y parameter', @ischar) addParameter(p, 'YLabelString', 'Y parameter', @ischar)
addParameter(p, 'TitleString', 'One-Parameter Scan', @ischar) addParameter(p, 'TitleString', 'One-Parameter Scan', @ischar)
@ -25,10 +22,7 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
YQuantity = p.Results.QuantityOfInterestArray; YQuantity = p.Results.QuantityOfInterestArray;
ErrorsForYQuantity = p.Results.ErrorsForYQuantity; ErrorsForYQuantity = p.Results.ErrorsForYQuantity;
ErrorsArray = p.Results.ErrorsArray; ErrorsArray = p.Results.ErrorsArray;
CIForYQuantity = p.Results.CIForYQuantity;
CIArray = p.Results.CIArray;
RescalingFactorForYQuantity = p.Results.RescalingFactorForYQuantity; RescalingFactorForYQuantity = p.Results.RescalingFactorForYQuantity;
RemoveOutliers = p.Results.RemoveOutliers;
YLabelString = p.Results.YLabelString; YLabelString = p.Results.YLabelString;
TitleString = p.Results.TitleString; TitleString = p.Results.TitleString;
@ -45,37 +39,16 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
screensize = get(0,'ScreenSize'); screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600]; f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
if RemoveOutliers
[YQuantity,TF] = rmoutliers(YQuantity);
XParameter = XParameter(~TF);
ErrorsArray = ErrorsArray(~TF);
ClippedCIArray(:,1) = CIArray(~TF,1);
ClippedCIArray(:,2) = CIArray(~TF,2);
CIArray = ClippedCIArray;
end
RescaledXParameter = XParameter .* RescalingFactorForXParameter; RescaledXParameter = XParameter .* RescalingFactorForXParameter;
RescaledYQuantity = YQuantity .* RescalingFactorForYQuantity; RescaledYQuantity = YQuantity .* RescalingFactorForYQuantity;
hold on
if ErrorsForYQuantity if ErrorsForYQuantity
RescaledErrorsArray = ErrorsArray .* RescalingFactorForYQuantity; RescaledErrorsArray = ErrorsArray .* RescalingFactorForYQuantity;
errorbar(RescaledXParameter, RescaledYQuantity, RescaledErrorsArray, 'o', 'Linewidth', 1.5, 'MarkerFaceColor', '#0071BB') errorbar(RescaledXParameter, RescaledYQuantity, RescaledErrorsArray)
else else
plot(RescaledXParameter, RescaledYQuantity, '--o', 'Linewidth', 1.5); plot(RescaledXParameter, RescaledYQuantity, '--o', 'Linewidth', 1.5);
end end
if CIForYQuantity
RescaledCIArray = CIArray .* RescalingFactorForYQuantity;
Plotting.plotConfidenceIntervalRegion(RescaledXParameter, RescaledCIArray(:,1), RescaledCIArray(:,2));
end
hold off
xlim([0 inf])
ylim([0 inf])
hXLabel = xlabel(XLabelString); hXLabel = xlabel(XLabelString);
hYLabel = ylabel(YLabelString); hYLabel = ylabel(YLabelString);
hTitle = sgtitle(TitleString); hTitle = sgtitle(TitleString);
@ -85,6 +58,5 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
set( hTitle , ... set( hTitle , ...
'FontSize' , 18 ); 'FontSize' , 18 );
grid on
Helper.bringFiguresWithTagInForeground(); Helper.bringFiguresWithTagInForeground();
end end

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

@ -24,18 +24,18 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
BlueDetuning; BlueDetuning;
BlueBeamRadius; BlueBeamRadius;
BlueBeamWaist; BlueBeamWaist;
BlueWaveNumber; BlueWaveVector;
BlueSaturationIntensity; BlueSaturationIntensity;
OrangePower; OrangePower;
OrangeDetuning; OrangeDetuning;
OrangeBeamRadius; OrangeBeamRadius;
OrangeBeamWaist; OrangeBeamWaist;
OrangeWaveNumber; OrangeWaveVector;
OrangeSaturationIntensity; OrangeSaturationIntensity;
CoolingBeamPower; CoolingBeamPower;
CoolingBeamWaveNumber; CoolingBeamWaveVector;
CoolingBeamLinewidth; CoolingBeamLinewidth;
CoolingBeamDetuning; CoolingBeamDetuning;
CoolingBeamRadius; CoolingBeamRadius;
@ -49,7 +49,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
SidebandBeamSaturationIntensity; SidebandBeamSaturationIntensity;
PushBeamPower; PushBeamPower;
PushBeamWaveNumber; PushBeamWaveVector;
PushBeamLinewidth; PushBeamLinewidth;
PushBeamDetuning; PushBeamDetuning;
PushBeamRadius; PushBeamRadius;
@ -71,16 +71,16 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
CaptureVelocity; CaptureVelocity;
VelocityCutoff; VelocityCutoff;
ClausingFactor; ClausingFactor;
ReducedClausingFactor; ThetaArray;
ReducedFlux; AngularDistribution;
TimeSpentInInteractionRegion; NormalizationConstantForAngularDistribution;
%Flags %Flags
SpontaneousEmission; SpontaneousEmission;
Sideband; Sideband;
ZeemanSlowerBeam; ZeemanSlowerBeam;
Gravity; Gravity;
BackgroundCollision; AtomicBeamCollision;
DebugMode; DebugMode;
DoSave; DoSave;
@ -124,7 +124,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
@islogical); @islogical);
addParameter(p, 'Gravity', false,... addParameter(p, 'Gravity', false,...
@islogical); @islogical);
addParameter(p, 'BackgroundCollision', false,... addParameter(p, 'AtomicBeamCollision', true,...
@islogical); @islogical);
addParameter(p, 'DebugMode', false,... addParameter(p, 'DebugMode', false,...
@islogical); @islogical);
@ -144,7 +144,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
s.Sideband = p.Results.Sideband; s.Sideband = p.Results.Sideband;
s.ZeemanSlowerBeam = p.Results.ZeemanSlowerBeam; s.ZeemanSlowerBeam = p.Results.ZeemanSlowerBeam;
s.Gravity = p.Results.Gravity; s.Gravity = p.Results.Gravity;
s.BackgroundCollision = p.Results.BackgroundCollision; s.AtomicBeamCollision = p.Results.AtomicBeamCollision;
s.DebugMode = p.Results.DebugMode; s.DebugMode = p.Results.DebugMode;
s.DoSave = p.Results.SaveData; s.DoSave = p.Results.SaveData;
@ -177,7 +177,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
ret = this.SimulationTime; ret = this.SimulationTime;
end end
function set.NumberOfAtoms(this, val) function set.NumberOfAtoms(this, val)
assert(val <= 20000, 'Not time efficient to compute for atom numbers larger than 20,000!'); assert(val <= 10000, 'Not time efficient to compute for atom numbers larger than 10,000!');
this.NumberOfAtoms = val; this.NumberOfAtoms = val;
end end
function ret = get.NumberOfAtoms(this) function ret = get.NumberOfAtoms(this)
@ -282,11 +282,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.BlueBeamWaist(this) function ret = get.BlueBeamWaist(this)
ret = this.BlueBeamWaist; ret = this.BlueBeamWaist;
end end
function set.BlueWaveNumber(this, val) function set.BlueWaveVector(this, val)
this.BlueWaveNumber = val; this.BlueWaveVector = val;
end end
function ret = get.BlueWaveNumber(this) function ret = get.BlueWaveVector(this)
ret = this.BlueWaveNumber; ret = this.BlueWaveVector;
end end
function set.BlueSaturationIntensity(this, val) function set.BlueSaturationIntensity(this, val)
this.BlueSaturationIntensity = val; this.BlueSaturationIntensity = val;
@ -319,11 +319,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.OrangeBeamWaist(this) function ret = get.OrangeBeamWaist(this)
ret = this.OrangeBeamWaist; ret = this.OrangeBeamWaist;
end end
function set.OrangeWaveNumber(this, val) function set.OrangeWaveVector(this, val)
this.OrangeWaveNumber = val; this.OrangeWaveVector = val;
end end
function ret = get.OrangeWaveNumber(this) function ret = get.OrangeWaveVector(this)
ret = this.OrangeWaveNumber; ret = this.OrangeWaveVector;
end end
function set.OrangeSaturationIntensity(this, val) function set.OrangeSaturationIntensity(this, val)
this.OrangeSaturationIntensity = val; this.OrangeSaturationIntensity = val;
@ -356,11 +356,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.CoolingBeamWaist(this) function ret = get.CoolingBeamWaist(this)
ret = this.CoolingBeamWaist; ret = this.CoolingBeamWaist;
end end
function set.CoolingBeamWaveNumber(this, val) function set.CoolingBeamWaveVector(this, val)
this.CoolingBeamWaveNumber = val; this.CoolingBeamWaveVector = val;
end end
function ret = get.CoolingBeamWaveNumber(this) function ret = get.CoolingBeamWaveVector(this)
ret = this.CoolingBeamWaveNumber; ret = this.CoolingBeamWaveVector;
end end
function set.CoolingBeamLinewidth(this, val) function set.CoolingBeamLinewidth(this, val)
this.CoolingBeamLinewidth = val; this.CoolingBeamLinewidth = val;
@ -430,11 +430,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.PushBeamWaist(this) function ret = get.PushBeamWaist(this)
ret = this.PushBeamWaist; ret = this.PushBeamWaist;
end end
function set.PushBeamWaveNumber(this, val) function set.PushBeamWaveVector(this, val)
this.PushBeamWaveNumber= val; this.PushBeamWaveVector = val;
end end
function ret = get.PushBeamWaveNumber(this) function ret = get.PushBeamWaveVector(this)
ret = this.PushBeamWaveNumber; ret = this.PushBeamWaveVector;
end end
function set.PushBeamLinewidth(this, val) function set.PushBeamLinewidth(this, val)
this.PushBeamLinewidth = val; this.PushBeamLinewidth = val;
@ -529,23 +529,31 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.ClausingFactor(this) function ret = get.ClausingFactor(this)
ret = this.ClausingFactor; ret = this.ClausingFactor;
end end
function set.ReducedClausingFactor(this,val) function set.AngularDistribution(this,val)
this.ReducedClausingFactor = val; this.AngularDistribution = val;
end end
function ret = get.ReducedClausingFactor(this) function ret = get.AngularDistribution(this)
ret = this.ReducedClausingFactor; ret = this.AngularDistribution;
end end
function set.ReducedFlux(this,val) function set.ThetaArray(this,val)
this.ReducedFlux = val; this.ThetaArray = val;
end end
function ret = get.ReducedFlux(this) function ret = get.ThetaArray(this)
ret = this.ReducedFlux; ret = this.ThetaArray;
end end
function set.TimeSpentInInteractionRegion(this,val) function set.NormalizationConstantForAngularDistribution(this,val)
this.TimeSpentInInteractionRegion = val; this.NormalizationConstantForAngularDistribution = val;
end end
function ret = get.TimeSpentInInteractionRegion(this) function ret = get.NormalizationConstantForAngularDistribution(this)
ret = this.TimeSpentInInteractionRegion; ret = this.NormalizationConstantForAngularDistribution;
end
function set.AtomicBeamCollision(this,val)
this.AtomicBeamCollision=val;
end
function ret=get.AtomicBeamCollision(this)
ret=this.AtomicBeamCollision;
end end
function set.DebugMode(this, val) function set.DebugMode(this, val)
@ -578,19 +586,19 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
methods methods
function ret = get.CoolingBeamSaturationParameter(this) function ret = get.CoolingBeamSaturationParameter(this)
ret = 0.1 * (4 * this.CoolingBeamPower) / (pi*this.CoolingBeamWaist^2 * this.CoolingBeamSaturationIntensity); % two beams are reflected ret = 4* this.CoolingBeamPower/(pi*this.CoolingBeamWaist^2)/this.CoolingBeamSaturationIntensity/10; % two beams are reflected
end end
function ret = get.SidebandSaturationParameter(this) function ret = get.SidebandSaturationParameter(this)
ret = 0.1 * (4 * this.SidebandPower) / (pi*this.SidebandBeamWaist^2 * this.SidebandBeamSaturationIntensity); ret = 4*this.SidebandPower/(pi*this.SidebandBeamWaist^2)/this.SidebandBeamSaturationIntensity/10;
end end
function ret = get.PushBeamSaturationParameter(this) function ret = get.PushBeamSaturationParameter(this)
ret = 0.1 * this.PushBeamPower/(pi * this.PushBeamWaist^2 * this.PushBeamSaturationIntensity); ret = this.PushBeamPower/(pi*this.PushBeamWaist^2)/this.PushBeamSaturationIntensity/10;
end end
function ret = get.ZeemanSlowerBeamSaturationParameter(this) function ret = get.ZeemanSlowerBeamSaturationParameter(this)
ret = 0.1 * this.ZeemanSlowerBeamPower / (pi * this.ZeemanSlowerBeamWaist^2 * this.ZeemanSlowerBeamSaturationIntensity); ret = this.ZeemanSlowerBeamPower/(pi*this.ZeemanSlowerBeamWaist^2)/this.ZeemanSlowerBeamSaturationIntensity/10;
end end
function ret = get.OvenTemperatureinKelvin(this) function ret = get.OvenTemperatureinKelvin(this)
@ -599,7 +607,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.AverageVelocity(this) function ret = get.AverageVelocity(this)
%See Background collision probability section in Barbiero %See Background collision probability section in Barbiero
ret = sqrt((8*Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (pi*Helper.PhysicsConstants.Dy164Mass)); ret = sqrt(((8*pi)/9) * ((Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)/Helper.PhysicsConstants.Dy164Mass));
end end
function ret = get.AtomicBeamDensity(this) function ret = get.AtomicBeamDensity(this)

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

@ -4,17 +4,17 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint); WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint);
Origin=[0,0,0]; Origin=[0,0,0];
SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist); SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationIntensity, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist);
DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveNumber; DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveVector;
Detuning = this.PushBeamDetuning - DopplerShift; Detuning = this.PushBeamDetuning - DopplerShift;
s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2)); a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveVector * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* ...
a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* (s_push/(1+s_push)); (SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2)));
if this.SpontaneousEmission if this.SpontaneousEmission
a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, Detuning, this.PushBeamLinewidth, this.PushBeamWaveNumber); a_scatter = this.accelerationDueToSpontaneousEmissionProcess(SaturationIntensity, SaturationIntensity, Detuning, this.PushBeamLinewidth, this.PushBeamWaveVector);
else else
a_scatter = [0,0,0]; a_scatter = [0,0,0];
end end

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

@ -1,12 +1,12 @@
function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationParameter, TotalSaturationParameter, Detuning, Linewidth, WaveNumber) function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationIntensity, TotalSaturationIntensity, Detuning, Linewidth, WaveVector)
Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1]; Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1];
Vector = Vector./norm(Vector); Vector = Vector./norm(Vector);
ScatteringRate = (0.5 * SaturationParameter * Linewidth) / (1 + TotalSaturationParameter + (4 * (Detuning / Linewidth)^2)); ScatteringRate = (0.5 * SaturationIntensity * Linewidth) / (1 + TotalSaturationIntensity + (4 * (Detuning/Linewidth)^2));
NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep); NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep);
if NumberOfScatteringEvents > 0 if NumberOfScatteringEvents > 0
ret = Vector.*((Helper.PhysicsConstants.PlanckConstantReduced * WaveNumber) / ... ret = Vector.*((Helper.PhysicsConstants.PlanckConstantReduced * WaveVector) / ...
(Helper.PhysicsConstants.Dy164Mass * this.TimeStep)).* sqrt(NumberOfScatteringEvents); (Helper.PhysicsConstants.Dy164Mass * this.TimeStep)).* sqrt(NumberOfScatteringEvents);
else else
ret = zeros(1,3); ret = zeros(1,3);

40
MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m
View File

@ -1,40 +0,0 @@
function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstimation(this, NumberOfLoadedAtoms)
n = this.NumberOfAtoms;
NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep);
Autocorrelation = autocorr(NumberOfLoadedAtoms,'NumLags', double(NumberOfTimeSteps - 1));
if Autocorrelation(1)~=0
CorrelationFactor = table(Helper.findAllZeroCrossings(linspace(1, double(NumberOfTimeSteps), double(NumberOfTimeSteps)), Autocorrelation)).Var1(1);
if ~isnan(CorrelationFactor)
SampleLength = floor(CorrelationFactor);
NumberOfBootsrapSamples = 1000;
MeanLoadingRatioInEachSample = zeros(1,NumberOfBootsrapSamples);
for SampleNumber = 1:NumberOfBootsrapSamples
BoostrapSample = datasample(NumberOfLoadedAtoms, SampleLength); % Sample with replacement
MeanLoadingRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
end
LoadingRate = mean(MeanLoadingRatioInEachSample) * this.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;
ts = tinv([0.025 0.975],NumberOfBootsrapSamples-1); % T-Score
ConfidenceInterval = LoadingRate + ts*StandardError; % 95% Confidence Intervals
else
LoadingRate = nan;
StandardError = nan;
ConfidenceInterval = [nan nan];
end
else
LoadingRate = nan;
StandardError = nan;
ConfidenceInterval = [nan nan];
end
end

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

@ -1,15 +1,14 @@
function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector) function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
switch this.SimulationMode VelocityVector = VelocityVector./norm(VelocityVector);
case "2D"
VelocityUnitVector = VelocityVector./norm(VelocityVector);
UpperLimit = 500; UpperLimit = 500;
LowerLimit = 0; LowerLimit = 0;
this.AtomicBeamCollision = false;
for Index = 1:500 for Index = 1:500
InitialVelocity = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; InitialVelocity = 0.5 * (UpperLimit + LowerLimit) * VelocityVector;
ParticleDynamicalQuantities = this.solver(PositionVector, InitialVelocity); [~, FinalDynamicalQuantities] = this.solver(PositionVector, InitialVelocity);
FinalPositionVector = ParticleDynamicalQuantities(end, 1:3); FinalPositionVector = FinalDynamicalQuantities(1:3);
if rssq(FinalPositionVector) <= this.OvenDistance if rssq(FinalPositionVector) <= this.OvenDistance
LowerLimit = 0.5 * (UpperLimit + LowerLimit); LowerLimit = 0.5 * (UpperLimit + LowerLimit);
else else
@ -17,12 +16,9 @@ function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
end end
if UpperLimit - LowerLimit < 1 if UpperLimit - LowerLimit < 1
ret = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; ret = InitialVelocity;
break; break;
end end
end end
case "3D"
% Development In progress
end
end end

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

@ -1,9 +0,0 @@
function ret = calculateClausingFactor(this)
ClausingFactorApproximation = (8 * this.NozzleRadius) / (3 * this.NozzleLength);
alpha = 0.5 - (3*this.Beta^2)^-1 * ((1 - (2*this.Beta^3) + ((2*this.Beta^2) - 1) * sqrt(1+this.Beta^2)) / (sqrt(1+this.Beta^2) - (this.Beta^2 * asinh((this.Beta^2)^-1))));
ClausingFactorAnalytic = 1 + (2/3) * (1 - (2 * alpha)) * (this.Beta - sqrt(1 - this.Beta^2)) + (2/3) * (1 + alpha) * this.Beta^(-2) * (1 - sqrt(1 + this.Beta^2));
ret = ClausingFactorAnalytic;
end

6
MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m
View File

@ -3,8 +3,6 @@ function ret = calculateFreeMolecularRegimeFlux(this)
%See Expected atomic flux section in Barbiero %See Expected atomic flux section in Barbiero
Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature
Dy164DensityinOven = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin); Dy164DensityinOven = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin);
ret = Helper.PhysicsConstants.Dy164IsotopicAbundance * 1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2; ClausingFactorApproximation = (8 * this.NozzleRadius) / (3 * this.NozzleLength);
% Needs to be multiplied with the "Clausing Factor" which here would be ret = Helper.PhysicsConstants.Dy164IsotopicAbundance * 1/4 * Dy164DensityinOven * this.AverageVelocity * ClausingFactorApproximation * pi * this.NozzleRadius.^2;
% the probability not for the full solid angle but the angle subtended
% by the aperture of the oven at its mouth.
end end

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

@ -1,38 +1,84 @@
function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities) function [LoadingRate, StandardError] = calculateLoadingRate(this, FinalDynamicalQuantities)
switch this.SimulationMode
case "2D"
n = this.NumberOfAtoms;
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 NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms);
for AtomIndex = 1:n AutocorrelationFunction = zeros(1, this.NumberOfAtoms);
TimeCounts(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
end for i = 1:this.NumberOfAtoms
this.TimeSpentInInteractionRegion = mean(TimeCounts); FinalPosition = FinalDynamicalQuantities(i,1:3);
for AtomIndex = 1:n DivergenceAngle = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2)));
CollisionEvents(AtomIndex) = this.computeCollisionProbability(); if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0)
if i == 1
NumberOfLoadedAtoms(1) = 1;
else
NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1) + 1;
end end
% Count the number of loaded atoms subject to conditions else
for TimeIndex = 1:NumberOfTimeSteps if i == 1
if TimeIndex ~= 1 NumberOfLoadedAtoms(1) = 0;
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1); else
end NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1);
for AtomIndex = 1:n
Position = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
Velocity = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
if this.exitCondition(Position, Velocity, CollisionEvents(AtomIndex))
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
end end
end end
end end
[LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(NumberOfLoadedAtoms); for i = 1:this.NumberOfAtoms-1
MeanTrappingEfficiency = 0;
MeanLoadingRateShifted = 0;
for j = 1:this.NumberOfAtoms-i
MeanTrappingEfficiency = MeanTrappingEfficiency + NumberOfLoadedAtoms(j)/j;
MeanLoadingRateShifted = MeanLoadingRateShifted + (NumberOfLoadedAtoms(i+j))/(i+j);
AutocorrelationFunction(i) = AutocorrelationFunction(i) + ((NumberOfLoadedAtoms(j)/j).*(NumberOfLoadedAtoms(i+j)/(i+j)));
end
AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanTrappingEfficiency * MeanLoadingRateShifted));
end
case "3D" if AutocorrelationFunction(1)~=0 % To handle cases where there is no atom loading
% Development In progress
AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1);
x = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms);
[FitObject,~] = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100);
CorrelationFactor = FitObject.n; % n is the autocorrelation factor
SumOfAllMeanTrappingEfficiencies = 0;
NumberOfBlocks = floor(this.NumberOfAtoms/(2*CorrelationFactor+1));
MeanTrappingEfficiencyInEachBlock = zeros(1,NumberOfBlocks);
BlockNumberLimit = min(NumberOfBlocks-1,5);
% Jackknife method is a systematic way of obtaining the “standard deviation” error of a
% set of stochastic measurements:
% 1. Calculate average (or some function f) from the full dataset.
% 2. Divide data into M blocks, with block length ≫ correlation factor . This is done in order to
% get rid of autocorrelations; if there are no correlations, block length can be 1.
% 3. For each m = 1 . . .M, take away block m and calculate the average using
% the data from all other blocks.
% 4. Estimate the error by calculating the deviation of the average in each block from average of the full dataset
%Jackknife estimate of the loading rate
for i = 1:NumberOfBlocks-BlockNumberLimit
MeanTrappingEfficiencyInEachBlock(i) = NumberOfLoadedAtoms(this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1))) / (this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1)));
SumOfAllMeanTrappingEfficiencies = SumOfAllMeanTrappingEfficiencies + MeanTrappingEfficiencyInEachBlock(i);
end end
MeanTrappingEfficiency = SumOfAllMeanTrappingEfficiencies / (NumberOfBlocks-BlockNumberLimit);
c = integral(@(velocity) sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ...
* velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* this.OvenTemperatureinKelvin))), 0, this.VelocityCutoff);
LoadingRate = MeanTrappingEfficiency * c * this.calculateFreeMolecularRegimeFlux();
% Jackknife estimate of the variance of the loading rate
Variance = 0;
for i = 1:NumberOfBlocks-BlockNumberLimit
Variance = Variance + (MeanTrappingEfficiencyInEachBlock(i) - MeanTrappingEfficiency)^2;
end
StandardError = sqrt((((NumberOfBlocks-BlockNumberLimit) - 1) / (NumberOfBlocks-BlockNumberLimit)) * Variance);
else
LoadingRate = 0;
StandardError = 0;
end
end end

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

@ -1,6 +1,12 @@
function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVector, WaveVectorOrigin, WaveVectorEndPoint, BeamRadius, BeamWaist) function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVector, WaveVectorOrigin, WaveVectorEndPoint, BeamRadius, BeamWaist)
WaveVector = WaveVectorEndPoint - WaveVectorOrigin; % Line
PositionVectorFromWaveVectorOrigin = PositionVector - WaveVectorOrigin; % Point = PositionVector
DistanceBetweenAtomAndLaserBeamAxis = Helper.calculateDistanceFromPointToLine(PositionVector, WaveVectorOrigin, WaveVectorEndPoint); %Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
%One side of parallelogram = PositionVectorFromWaveVectorOrigin
%Base = Wavevector
%Area = One side of parallelogram X Base
DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector);
if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius
ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2); ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2);

17
MOT Capture Process Simulation/@MOTSimulator/calculateReducedClausingFactor.m
View File

@ -1,17 +0,0 @@
function [ReducedClausingFactor, NormalizationConstantForAngularDistribution] = calculateReducedClausingFactor(this)
ThetaArray = linspace(0.0001, pi/2, 1000);
AngularDistribution = zeros(1,length(ThetaArray));
parfor k = 1:length(ThetaArray)
AngularDistribution(k) = this.angularDistributionFunction(ThetaArray(k));
end
NormalizationConstantForAngularDistribution = max(2 * pi .* sin(ThetaArray) .* AngularDistribution);
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
ReducedClausingFactor = ReducedClausingFactor + (2 * pi * sin(ThetaArray(p)) * AngularDistribution(p) * (ThetaArray(2)-ThetaArray(1)));
end
end
ReducedClausingFactor = ReducedClausingFactor / pi;
end

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

@ -28,9 +28,9 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
B = sign(LocalMagneticField(1:3) * WaveVectorEndPoint(i,1:3)') * LocalMagneticField(4); B = sign(LocalMagneticField(1:3) * WaveVectorEndPoint(i,1:3)') * LocalMagneticField(4);
ZeemanShift = this.LandegFactor * this.MagneticSubLevel * (Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced) * B; ZeemanShift = this.LandegFactor * this.MagneticSubLevel * Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced * B;
DopplerShift = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveNumber; DopplerShift = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveVector;
Delta_Cooling(i*2-1) = this.CoolingBeamDetuning + DopplerShift + ZeemanShift * Sigma(i); 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) = this.CoolingBeamDetuning - DopplerShift - ZeemanShift * Sigma(i);
@ -44,11 +44,11 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
SaturationParameter = [0,0,0,0,0,0,0,0]; SaturationParameter = [0,0,0,0,0,0,0,0];
for i = 1:2 for i = 1:2
SaturationParameter(2*i-1) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2); SaturationParameter(2*i-1) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4*Delta_Cooling(2*i-1)^2 / this.CoolingBeamLinewidth^2);
SaturationParameter(2*i) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Cooling(2*i) /this.CoolingBeamLinewidth)^2); SaturationParameter(2*i) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4*Delta_Cooling(2*i)^2 / this.CoolingBeamLinewidth^2);
if this.Sideband if this.Sideband
SaturationParameter(2*i-1+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2); SaturationParameter(2*i-1+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(1)) /(1 + 4*Delta_Sideband(2*i-1)^2/ this.CoolingBeamLinewidth^2);
SaturationParameter(2*i+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2)) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2); SaturationParameter(2*i+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2)) /(1 + 4*Delta_Sideband(2*i)^2 / this.CoolingBeamLinewidth^2);
end end
end end
@ -56,28 +56,28 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
for i = 1:2 for i = 1:2
a_1 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... a_1 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
(SaturationParameter(2*i-1)/(1 + TotalSaturationParameter)); (SaturationParameter(2*i-1)/(1 + TotalSaturationParameter));
a_2 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... a_2 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
(SaturationParameter(2*i) / (1 + TotalSaturationParameter)); (SaturationParameter(2*i) / (1 + TotalSaturationParameter));
if this.SpontaneousEmission if this.SpontaneousEmission
a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ... a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveVector) + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, Delta_Cooling(2*i), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber); this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, Delta_Cooling(2*i) , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector);
else else
a_scattering = [0,0,0]; a_scattering = [0,0,0];
end end
if this.Sideband if this.Sideband
a_1 = a_1 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... a_1 = a_1 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
(SaturationParameter(2*i-1+4)/(1 + TotalSaturationParameter)); (SaturationParameter(2*i-1+4)/(1 + TotalSaturationParameter));
a_2 = a_2 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... a_2 = a_2 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
(SaturationParameter(2*i+4)/(1 + TotalSaturationParameter)); (SaturationParameter(2*i+4)/(1 + TotalSaturationParameter));
if this.SpontaneousEmission if this.SpontaneousEmission
a_scattering = a_scattering + ... a_scattering = a_scattering + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ... this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveVector) + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4), TotalSaturationParameter, Delta_Cooling(2*i), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber); this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4), TotalSaturationParameter, Delta_Cooling(2*i) , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector);
else else
a_scattering = [0,0,0]; a_scattering = [0,0,0];
end end

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

@ -1,9 +0,0 @@
function ret = computeCollisionProbability(this)
collision = rand(1);
CollisionProbability = 1 - exp(-this.TimeSpentInInteractionRegion/this.CollisionTime);
if this.BackgroundCollision && collision <= CollisionProbability
ret = true;
else
ret = false;
end
end

20
MOT Capture Process Simulation/@MOTSimulator/computeTimeSpentInInteractionRegion.m
View File

@ -1,20 +0,0 @@
function T = computeTimeSpentInInteractionRegion(this, r)
% INPUT:
% r : N x 3 array. N is the number of time steps
% OUTPUT
% T : gives the distribution of time spent in the interaction region
% USAGE:
% T = this.computeTimeSpentInInteractionRegion(r)
T = 0;
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
A = 1;
else
A = 0;
end
T = T + A * this.TimeStep;
end
end

8
MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m
View File

@ -1,4 +1,4 @@
function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOneParameterScan(this, ParameterName, ParameterArray, varargin) function [LoadingRateArray, StandardErrorArray] = doOneParameterScan(this, ParameterName, ParameterArray, varargin)
p = inputParser; p = inputParser;
p.KeepUnmatched = true; p.KeepUnmatched = true;
@ -30,7 +30,6 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOne
NumberOfPointsForParam = length(ParameterArray); NumberOfPointsForParam = length(ParameterArray);
LoadingRateArray = zeros(1,NumberOfPointsForParam); LoadingRateArray = zeros(1,NumberOfPointsForParam);
StandardErrorArray = zeros(1,NumberOfPointsForParam); StandardErrorArray = zeros(1,NumberOfPointsForParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2);
for i=1:NumberOfPointsForParam for i=1:NumberOfPointsForParam
eval(sprintf('OptionsStruct.%s = %d;', ParameterName, ParameterArray(i))); eval(sprintf('OptionsStruct.%s = %d;', ParameterName, ParameterArray(i)));
@ -45,18 +44,15 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOne
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
this.setInitialConditions(options{:}); this.setInitialConditions(options{:});
tic tic
[LoadingRate, StandardError, ConfidenceInterval] = this.runSimulation(); [LoadingRate, StandardError] = this.runSimulation();
LoadingRateArray(i) = LoadingRate; LoadingRateArray(i) = LoadingRate;
StandardErrorArray(i) = StandardError; StandardErrorArray(i) = StandardError;
ConfidenceIntervalArray(i,1) = ConfidenceInterval(1);
ConfidenceIntervalArray(i,2) = ConfidenceInterval(2);
end end
if this.DoSave if this.DoSave
LoadingRate = struct; LoadingRate = struct;
LoadingRate.Values = LoadingRateArray; LoadingRate.Values = LoadingRateArray;
LoadingRate.Errors = StandardErrorArray; LoadingRate.Errors = StandardErrorArray;
LoadingRate.CI = ConfidenceIntervalArray;
this.Results = LoadingRate; this.Results = LoadingRate;
SaveFolder = [this.SaveDirectory filesep 'Results']; SaveFolder = [this.SaveDirectory filesep 'Results'];
Filename = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')]; Filename = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')];

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

@ -1,4 +1,4 @@
function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ... function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ...
SecondParameterName, SecondParameterArray, varargin) SecondParameterName, SecondParameterArray, varargin)
p = inputParser; p = inputParser;
@ -38,7 +38,6 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwo
NumberOfPointsForSecondParam = length(SecondParameterArray); NumberOfPointsForSecondParam = length(SecondParameterArray);
LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
for i=1:NumberOfPointsForFirstParam for i=1:NumberOfPointsForFirstParam
eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i))); eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i)));
@ -59,11 +58,9 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwo
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
this.setInitialConditions(options{:}); this.setInitialConditions(options{:});
tic tic
[LoadingRate, StandardError, ConfidenceInterval] = this.runSimulation(); [LoadingRate, StandardError] = this.runSimulation();
LoadingRateArray(i,j) = LoadingRate; LoadingRateArray(i,j) = LoadingRate;
StandardErrorArray(i,j) = StandardError; StandardErrorArray(i,j) = StandardError;
ConfidenceIntervalArray(i, j, 1) = ConfidenceInterval(1);
ConfidenceIntervalArray(i, j, 2) = ConfidenceInterval(2);
end end
end end
@ -71,7 +68,6 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwo
LoadingRate = struct; LoadingRate = struct;
LoadingRate.Values = LoadingRateArray; LoadingRate.Values = LoadingRateArray;
LoadingRate.Errors = StandardErrorArray; LoadingRate.Errors = StandardErrorArray;
LoadingRate.CI = ConfidenceIntervalArray;
this.Results = LoadingRate; this.Results = LoadingRate;
SaveFolder = [this.SaveDirectory filesep 'Results']; SaveFolder = [this.SaveDirectory filesep 'Results'];
Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')]; Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];

11
MOT Capture Process Simulation/@MOTSimulator/exitCondition.m
View File

@ -1,11 +0,0 @@
function ret = exitCondition(this, 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
ret = true;
else
ret = false;
end
end

5
MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m
View File

@ -1,12 +1,7 @@
function ret = initialPositionSampling(this) function ret = initialPositionSampling(this)
switch this.SimulationMode
case "2D"
n = this.NumberOfAtoms; n = this.NumberOfAtoms;
phi = 2 * pi * rand(n,1); phi = 2 * pi * rand(n,1);
rho = this.Beta * 0.5 * this.NozzleLength * sqrt(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)]; ret = [-this.OvenDistance * ones(n,1), rho.*cos(phi), rho.*sin(phi)];
case "3D"
% Development In progress
end
end end

42
MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m
View File

@ -1,25 +1,6 @@
function ret = initialVelocitySampling(this) function ret = initialVelocitySampling(this)
switch this.SimulationMode
case "2D"
n = this.NumberOfAtoms; 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)/integral(VelocityDistribution,0,Inf);
this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
ret = zeros(n,3); ret = zeros(n,3);
SampledVelocityMagnitude = zeros(n,1); SampledVelocityMagnitude = zeros(n,1);
SampledPolarAngle = zeros(n,1); SampledPolarAngle = zeros(n,1);
@ -32,25 +13,27 @@ function ret = initialVelocitySampling(this)
else else
MaximumVelocityAllowed = MostProbableVelocity; MaximumVelocityAllowed = MostProbableVelocity;
end end
NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed); ProbabilityOfMaximumVelocityAllowed = this.velocityDistributionFunction(MaximumVelocityAllowed);
ProbabilityOfMaximumDivergenceAngleAllowed = 0.98 * this.NormalizationConstantForAngularDistribution * max(this.AngularDistribution .* sin(this.ThetaArray));
parfor i = 1:n parfor i = 1:n
% Rejection Sampling of speed % Rejection Sampling of speed
y = rand(1); y = ProbabilityOfMaximumVelocityAllowed * rand(1);
x = this.VelocityCutoff * 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); while y > this.velocityDistributionFunction(x) %As long as this loop condition is satisfied, reject the corresponding x value
y = ProbabilityOfMaximumVelocityAllowed * rand(1);
x = this.VelocityCutoff * rand(1); x = this.VelocityCutoff * rand(1);
end end
SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample
% Rejection Sampling of polar angle % Rejection Sampling of polar angle
w = rand(1); z = this.MOTExitDivergence * rand(1);
z = this.NozzleExitDivergence * rand(1); w = ProbabilityOfMaximumDivergenceAngleAllowed * 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 while w > (this.NormalizationConstantForAngularDistribution * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value
w = rand(1); z = this.MOTExitDivergence * rand(1);
z = this.NozzleExitDivergence * rand(1); w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
end end
SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
@ -60,8 +43,5 @@ function ret = initialVelocitySampling(this)
ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ... ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ...
SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))]; SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))];
end end
case "3D"
% Development In progress
end
end end

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

@ -3,9 +3,8 @@ function reinitializeSimulator(this)
pc = Helper.PhysicsConstants; pc = Helper.PhysicsConstants;
%% SIMULATION PARAMETERS %% SIMULATION PARAMETERS
this.NozzleLength = 60e-3; this.NozzleLength = 60e-3;
this.NozzleRadius = 2.60e-3; this.NozzleRadius = 2.50e-3;
this.Beta = 2 * (this.NozzleRadius/this.NozzleLength); this.Beta = 2 * (this.NozzleRadius/this.NozzleLength);
this.ClausingFactor = this.calculateClausingFactor();
this.ApertureCut = max(2.5e-3,this.NozzleRadius); this.ApertureCut = max(2.5e-3,this.NozzleRadius);
this.OvenDistance = (25+12.5)*1e-3 + (this.NozzleRadius + this.ApertureCut) / tan(15/360 * 2 * pi); this.OvenDistance = (25+12.5)*1e-3 + (this.NozzleRadius + this.ApertureCut) / tan(15/360 * 2 * pi);
% Distance between the nozzle and the 2-D MOT chamber center % Distance between the nozzle and the 2-D MOT chamber center
@ -13,25 +12,20 @@ function reinitializeSimulator(this)
% 12.5 is the radius of the oven % 12.5 is the radius of the oven
% 15 eg is the angle between the 2-D MOT chamber center and the nozzle % 15 eg is the angle between the 2-D MOT chamber center and the nozzle
this.OvenTemperature = 1000; % Temperature in Celsius this.OvenTemperature = 1000; % Temperature in Celsius
this.MOTDistance = 0.32; % Distance between the 2-D MOT the 3-D MOT this.MOTDistance = 320e-3; % Distance between the 2-D MOT the 3-D MOT
this.BlueWaveNumber = 2*pi/pc.BlueWavelength; this.BlueWaveVector = 2*pi/pc.BlueWavelength;
this.BlueSaturationIntensity = 0.1 * (2 * pi^2 / 3) * ((pc.PlanckConstantReduced * pc.SpeedOfLight * pc.BlueLinewidth) / (pc.BlueWavelength)^3); this.BlueSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.BlueLinewidth/3/(pc.BlueWavelength)^3/10;
this.OrangeWaveNumber = 2*pi/pc.OrangeWavelength; this.OrangeWaveVector = 2*pi/pc.OrangeWavelength;
this.OrangeSaturationIntensity = 0.1 * (2 * pi^2 / 3) * ((pc.PlanckConstantReduced * pc.SpeedOfLight * pc.OrangeLinewidth) / (pc.OrangeWavelength)^3); this.OrangeSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.OrangeLinewidth/3/(pc.OrangeWavelength)^3/10;
this.BlueBeamRadius = min(0.035/2,sqrt(2)/2*this.OvenDistance); % Diameter of CF40 flange = 0.035 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_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 Theta_Aperture = 15/360*2*pi; % The limitation angle of the second aperture in the oven
this.NozzleExitDivergence = min(Theta_Nozzle,Theta_Aperture); this.NozzleExitDivergence = min(Theta_Nozzle,Theta_Aperture);
this.MOTExitDivergence = 16e-3; % The limitation angle between 2D-MOT and 3D-MOT this.MOTExitDivergence = 0.016; % The limitation angle between 2D-MOT and 3D-MOT
this.TimeSpentInInteractionRegion = this.SimulationTime; this.TotalPower = 0.4;
this.TotalPower = 0.8; this.OrangeBeamRadius = 1.2e-03;
this.OrangeBeamRadius = 1.2e-3; this.PushBeamRadius = 0.000;
this.PushBeamRadius = 1.2e-3;
this.PushBeamDistance = 0.32; 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; this.DistanceBetweenPushBeamAnd3DMOTCenter = 0;
this.ZeemanSlowerBeamRadius = 1;
end end

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

@ -1,6 +1,36 @@
function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this) function [LoadingRate, StandardError] = runSimulation(this)
n = this.NumberOfAtoms;
%% Calculate Background Collision Time --> Calculate Capture velocity --> Introduce velocity cutoff --> Calculate capture fraction
%% - Sampling for initial positions and velocities this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); % Make 5% larger than the numerically obtained CV
this.VelocityCutoff = this.CaptureVelocity(1); %Should be the magnitude of the 3-D velocity vector but since here the obtained capture
%velocity is only along the x-axis, we take the first term which is the x-component of the velocity.
%% Calculate the Clausing Factor
% Compute the angular distribution of the atomic beam
ThetaArray = linspace(0.0001, pi/2, 1000);
AngularDistribution = zeros(1,length(ThetaArray));
parfor i = 1:length(ThetaArray)
AngularDistribution(i) = this.angularDistributionFunction(ThetaArray(i));
end
% Numerically integrate the angular distribution over the full solid angle
NormalizationConstant = 0;
for j = 1:length(ThetaArray)
if ThetaArray(j) <= this.NozzleExitDivergence
NormalizationConstant = NormalizationConstant + (2 * pi * sin(ThetaArray(j)) * AngularDistribution(j) * (ThetaArray(2)-ThetaArray(1)));
end
end
this.ClausingFactor = (1/pi) * NormalizationConstant; %The complete intergration will give pi * ClausingFactor.
%Therefore, the Clausing Factor needs to be extracted from the
%result of the above integration by dividing out pi
this.ThetaArray = ThetaArray;
this.AngularDistribution = AngularDistribution;
this.NormalizationConstantForAngularDistribution = NormalizationConstant;
%%
% - sampling the position distribution % - sampling the position distribution
this.InitialPositions = this.initialPositionSampling(); this.InitialPositions = this.initialPositionSampling();
% - sampling the velocity distribution % - sampling the velocity distribution
@ -8,18 +38,18 @@ function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this)
%% Solve ODE %% Solve ODE
progressbar = Helper.parforNotifications(); progressbar = Helper.parforNotifications();
progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']); progressbar.PB_start(n,'Message',['Simulating capture process for ' num2str(n,'%.0f') ' atoms:']);
% calculate the final position of the atoms % calculate the final position of the atoms
ParticleDynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6); FinalDynamicalQuantities = zeros(n,9);
Positions = this.InitialPositions; Positions = this.InitialPositions;
Velocities = this.InitialVelocities; Velocities = this.InitialVelocities;
parfor Index = 1:this.NumberOfAtoms parfor Index = 1:n
ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:)); ret = this.solver(Positions(Index,:), Velocities(Index,:));
FinalDynamicalQuantities(Index,:) = ret(end,:);
progressbar.PB_iterate(); progressbar.PB_iterate();
end end
clear Index clear Index
%% Calculate the Loading Rate %% Calculate the Loading Rate
[LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ParticleDynamicalQuantities); [LoadingRate, StandardError] = this.calculateLoadingRate(FinalDynamicalQuantities);
end end

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

@ -4,27 +4,27 @@ function setInitialConditions(this, varargin)
p.KeepUnmatched = true; p.KeepUnmatched = true;
addParameter(p, 'NumberOfAtoms', 5000,... addParameter(p, 'NumberOfAtoms', 5000,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'BluePower', this.TotalPower,... addParameter(p, 'BluePower', 200e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'BlueDetuning', -1.39*Helper.PhysicsConstants.BlueLinewidth,... addParameter(p, 'BlueDetuning', -1.5*Helper.PhysicsConstants.BlueLinewidth,...
@(x) assert(isnumeric(x) && isscalar(x))); @(x) assert(isnumeric(x) && isscalar(x)));
addParameter(p, 'BlueBeamWaist', 16.6667e-3,... addParameter(p, 'BlueBeamWaist', 10e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'SidebandPower', 0.5*this.TotalPower,... addParameter(p, 'SidebandPower', 200e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'SidebandDetuning', -4.5*Helper.PhysicsConstants.BlueLinewidth,... addParameter(p, 'SidebandDetuning', 0,...
@(x) assert(isnumeric(x) && isscalar(x))); @(x) assert(isnumeric(x) && isscalar(x)));
addParameter(p, 'SidebandBeamWaist', 15e-3,... addParameter(p, 'SidebandBeamWaist', 12e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'PushBeamPower', 100e-3,... addParameter(p, 'PushBeamPower', 10e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'PushBeamDetuning', -5*Helper.PhysicsConstants.OrangeLinewidth,... addParameter(p, 'PushBeamDetuning', 0,...
@(x) assert(isnumeric(x) && isscalar(x))); @(x) assert(isnumeric(x) && isscalar(x)));
addParameter(p, 'PushBeamWaist', 0.81e-3,... addParameter(p, 'PushBeamWaist', 0.81e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'OrangePower', 70e-3,... addParameter(p, 'OrangePower', 70e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'OrangeDetuning', -1*Helper.PhysicsConstants.OrangeLinewidth,... addParameter(p, 'OrangeDetuning', 12e-3,...
@(x) assert(isnumeric(x) && isscalar(x))); @(x) assert(isnumeric(x) && isscalar(x)));
addParameter(p, 'OrangeBeamWaist', 12e-3,... addParameter(p, 'OrangeBeamWaist', 12e-3,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
@ -56,19 +56,24 @@ function setInitialConditions(this, varargin)
this.ZeemanSlowerBeamDetuning = p.Results.ZeemanSlowerBeamDetuning; this.ZeemanSlowerBeamDetuning = p.Results.ZeemanSlowerBeamDetuning;
this.ZeemanSlowerBeamWaist = p.Results.ZeemanSlowerBeamWaist; this.ZeemanSlowerBeamWaist = p.Results.ZeemanSlowerBeamWaist;
this.MagneticGradient = p.Results.MagneticGradient; this.MagneticGradient = p.Results.MagneticGradient;
%% Set general parameters according to simulation mode %% Set general parameters according to simulation mode
switch this.SimulationMode switch this.SimulationMode
case "2D" case "2D"
this.CoolingBeamPower = this.BluePower; this.CoolingBeamPower = this.BluePower;
this.CoolingBeamWaist = this.BlueBeamWaist; this.CoolingBeamWaist = this.BlueBeamWaist;
this.CoolingBeamLinewidth = Helper.PhysicsConstants.BlueLinewidth; this.CoolingBeamLinewidth = Helper.PhysicsConstants.BlueLinewidth;
this.CoolingBeamWaveNumber = this.BlueWaveNumber; this.CoolingBeamWaveVector = this.BlueWaveVector;
this.CoolingBeamDetuning = this.BlueDetuning; this.CoolingBeamDetuning = this.BlueDetuning;
this.CoolingBeamRadius = this.BlueBeamRadius; this.CoolingBeamRadius = this.BlueBeamRadius;
this.CoolingBeamWaist = this.BlueBeamWaist; this.CoolingBeamWaist = this.BlueBeamWaist;
this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity; this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity;
this.SidebandBeamRadius = this.BlueBeamRadius; this.SidebandBeamRadius = this.BlueBeamRadius;
this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity; this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity;
this.PushBeamLinewidth = Helper.PhysicsConstants.OrangeLinewidth;
this.PushBeamWaveVector = this.OrangeWaveVector;
this.PushBeamDetuning = this.OrangeDetuning;
this.PushBeamSaturationIntensity = this.OrangeSaturationIntensity;
this.LandegFactor = Helper.PhysicsConstants.BlueLandegFactor; this.LandegFactor = Helper.PhysicsConstants.BlueLandegFactor;
this.MagneticSubLevel = 1; this.MagneticSubLevel = 1;
case "3D" case "3D"

25
MOT Capture Process Simulation/@MOTSimulator/solver.m
View File

@ -1,16 +1,22 @@
function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelocity) function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity)
if this.Gravity if this.Gravity
g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration]; g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration];
else else
g = 0; g = 0;
end end
ParticleDynamicalQuantities = zeros(int64(this.SimulationTime/this.TimeStep),6); % Probability of Background Collisions
collision = rand(1);
CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime);
if collision >= CollisionProbability || this.AtomicBeamCollision == false %flag for skipping the background collision
ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9);
for i=1:int64(this.SimulationTime/this.TimeStep) for i=1:int64(this.SimulationTime/this.TimeStep)
ParticleDynamicalQuantities(i,1:3) = InitialPosition; ParticleTrajectory(i,1:3) = InitialPosition;
ParticleDynamicalQuantities(i,4:6) = InitialVelocity; ParticleTrajectory(i,4:6) = InitialVelocity;
rt = InitialPosition; rt = InitialPosition;
vt = InitialVelocity; vt = InitialVelocity;
@ -35,7 +41,14 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelo
InitialPosition = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6; InitialPosition = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6;
InitialVelocity = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6; InitialVelocity = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6;
%Acceleration = (ga1+2*(ga2+ga3)+ ga4)./6; ParticleTrajectory(i,7:9) = (ga1+2*(ga2+ga3)+ ga4)./6;
end
FinalDynamicalQuantities = ParticleTrajectory(end,:);
else
ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9);
FinalDynamicalQuantities = -500*ones(1,9); % -500 is a flag for giving up this trajectory
end end
end end

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

@ -1,5 +1,5 @@
function ret = velocityDistributionFunction(this, velocity) function ret = velocityDistributionFunction(this, velocity)
ret = sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin))^3 ... ret = sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ...
* velocity^3 * exp(-velocity^2.*(Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ... * velocity^3 * exp(-velocity^2.*(Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* this.OvenTemperatureinKelvin))); * this.OvenTemperatureinKelvin)));
end end

58
MOT Capture Process Simulation/test_MOTSimulator.m
View File

@ -7,10 +7,12 @@ OptionsStruct.SimulationMode = '2D';
OptionsStruct.TimeStep = 50e-06; % in s OptionsStruct.TimeStep = 50e-06; % in s
OptionsStruct.SimulationTime = 4e-03; % in s OptionsStruct.SimulationTime = 4e-03; % in s
OptionsStruct.SpontaneousEmission = true; OptionsStruct.SpontaneousEmission = true;
OptionsStruct.Sideband = false; OptionsStruct.Sideband = true;
OptionsStruct.ZeemanSlowerBeam = false;
OptionsStruct.Gravity = true; OptionsStruct.Gravity = true;
OptionsStruct.BackgroundCollision = true; OptionsStruct.AtomicBeamCollision = true;
OptionsStruct.SaveData = false; OptionsStruct.DebugMode = false;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation'; OptionsStruct.SaveDirectory = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation';
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
@ -40,18 +42,9 @@ clear OptionsStruct
[LoadingRate, ~] = Simulator.runSimulation(); [LoadingRate, ~] = Simulator.runSimulation();
%% - Plot initial distribution %% - Plot initial distribution
Simulator.setInitialConditions();
% - sampling the position distribution
Simulator.InitialPositions = Simulator.initialPositionSampling();
% - sampling the velocity distribution
Simulator.InitialVelocities = Simulator.initialVelocitySampling();
NumberOfBins = 100; NumberOfBins = 100;
Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins); Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins);
%% - Plot distributions of magnitude and direction of initial velocities
NumberOfBins = 50;
Plotting.plotInitialVeloctiySamplingVsAngle(Simulator, NumberOfBins)
%% - Plot Magnetic Field %% - Plot Magnetic Field
XAxisRange = [-5 5]; XAxisRange = [-5 5];
YAxisRange = [-5 5]; YAxisRange = [-5 5];
@ -63,7 +56,7 @@ TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius
Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius) Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
%% - Plot the Free Molecular Flux for different temperatures %% - Plot the Free Molecular Flux for different temperatures
Temperature = [950, 1000, 1050]; % Temperature Temperature = [900, 950, 1000, 1050, 1100]; % Temperature'
Plotting.plotFreeMolecularFluxVsTemp(Simulator,Temperature) Plotting.plotFreeMolecularFluxVsTemp(Simulator,Temperature)
%% - Plot Angular Distribution for different Beta %% - Plot Angular Distribution for different Beta
@ -71,12 +64,18 @@ Beta = [0.5, 0.1 , 0.05, 0.02, 0.01]; %Beta = 2 * radius / length of the tube
Plotting.plotAngularDistributionForDifferentBeta(Simulator, Beta) Plotting.plotAngularDistributionForDifferentBeta(Simulator, Beta)
%% - Plot Capture Velocity %% - Plot Capture Velocity
Simulator.setInitialConditions(); Simulator.setInitialConditions();
Simulator.SimulationTime = 15*10^(-4);
Simulator.TimeStep = 1.5*10^(-6);
Simulator.MOTExitDivergence = 15/360*2*pi;
Simulator.MagneticGradient = 0.4;
Plotting.plotCaptureVelocityVsAngle(Simulator); Plotting.plotCaptureVelocityVsAngle(Simulator);
%% - Plot Phase Space with Acceleration Field %% - Plot Phase Space with Acceleration Field
Simulator.NumberOfAtoms = 200; Simulator.NumberOfAtoms = 100;
MaximumVelocity = 150; MaximumVelocity = 150;
NumberOfBins = 200; %Along each axis NumberOfBins = 200; %Along each axis
IncidentAtomDirection = 0*2*pi/360; IncidentAtomDirection = 0*2*pi/360;
@ -87,7 +86,7 @@ Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberO
% ONE-PARAMETER SCAN % ONE-PARAMETER SCAN
NumberOfPointsForParam = 10; %iterations of the simulation NumberOfPointsForParam = 20; %iterations of the simulation
% Scan Cooling Beam Power % Scan Cooling Beam Power
PowerArray = linspace(0.1, 1.0, NumberOfPointsForParam) * Simulator.TotalPower; PowerArray = linspace(0.1, 1.0, NumberOfPointsForParam) * Simulator.TotalPower;
% Scan Cooling Beam Detuning % Scan Cooling Beam Detuning
@ -96,12 +95,11 @@ PowerArray = linspace(0.1, 1.0, NumberOfPointsForParam) * Simulator.TotalPow
OptionsStruct = struct; OptionsStruct = struct;
OptionsStruct.ChangeInitialConditions = true; OptionsStruct.ChangeInitialConditions = true;
OptionsStruct.ParameterNameArray = {'NumberOfAtoms'}; OptionsStruct.ParameterNameArray = {'NumberOfAtoms'};
OptionsStruct.ParameterValueArray = {5000}; OptionsStruct.ParameterValueArray = {10000};
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
tStart = tic; tStart = tic;
[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:}); [LoadingRateArray, ~] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
tEnd = toc(tStart); tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
@ -109,19 +107,15 @@ clear OptionsStruct
% - Plot results % - Plot results
ParameterArray = PowerArray; ParameterArray = PowerArray .* Simulator.TotalPower;
QuantityOfInterestArray = LoadingRateArray; QuantityOfInterestArray = LoadingRateArray;
OptionsStruct = struct; OptionsStruct = struct;
OptionsStruct.RescalingFactorForParameter = 1000; OptionsStruct.RescalingFactorForParameter = 1000;
OptionsStruct.XLabelString = 'Cooling Beam Power (mW)'; OptionsStruct.XLabelString = 'Cooling Beam Power (mW)';
OptionsStruct.RescalingFactorForYQuantity = 1e-10; OptionsStruct.RescalingFactorForYQuantity = 1e-16;
OptionsStruct.ErrorsForYQuantity = true; OptionsStruct.ErrorsForYQuantity = false;
OptionsStruct.ErrorsArray = StandardErrorArray; OptionsStruct.YLabelString = 'Loading rate (x 10^{16} atoms/s)';
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)', Simulator.MagneticGradient * 100); OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
@ -132,11 +126,11 @@ clear OptionsStruct
%% TWO-PARAMETER SCAN %% TWO-PARAMETER SCAN
NumberOfPointsForParam = 10; %iterations of the simulation NumberOfPointsForParam = 20; %iterations of the simulation
NumberOfPointsForSecondParam = 10; NumberOfPointsForSecondParam = 10;
% Scan Sideband Detuning and Power Ratio % Scan Sideband Detuning and Power Ratio
DetuningArray = linspace(-0.5,-10, NumberOfPointsForParam) * Helper.PhysicsConstants.BlueLinewidth; DetuningArray = linspace(-0.5,-6, NumberOfPointsForParam) * Helper.PhysicsConstants.BlueLinewidth;
SidebandPowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam) * Simulator.TotalPower; SidebandPowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam) * Simulator.TotalPower;
BluePowerArray = Simulator.TotalPower - SidebandPowerArray; BluePowerArray = Simulator.TotalPower - SidebandPowerArray;
@ -147,11 +141,11 @@ OptionsStruct.RelatedParameterName = 'BluePower';
OptionsStruct.RelatedParameterArray = BluePowerArray; OptionsStruct.RelatedParameterArray = BluePowerArray;
OptionsStruct.ChangeInitialConditions = true; OptionsStruct.ChangeInitialConditions = true;
OptionsStruct.ParameterNameArray = {'NumberOfAtoms'}; OptionsStruct.ParameterNameArray = {'NumberOfAtoms'};
OptionsStruct.ParameterValueArray = {5000}; OptionsStruct.ParameterValueArray = {10000};
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
tStart = tic; tStart = tic;
[LoadingRateArray, StandardErrorArray, ConfidenceInterval] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:}); [LoadingRateArray, StandardErrorArray] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
tEnd = toc(tStart); tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
@ -168,8 +162,8 @@ OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.B
OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)'; OptionsStruct.XLabelString = 'Sideband Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000; OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Sideband Power (mW)'; OptionsStruct.YLabelString = 'Sideband Power (mW)';
OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-10; OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-16;
OptionsStruct.ZLabelString = 'Loading rate (x 10^{10} atoms/s)'; OptionsStruct.ZLabelString = 'Loading rate (x 10^{16} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100); OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);