Compare commits
No commits in common. "c2e6420243e31a3385f4b9a311aa7ca60f0eae9c" and "d0b4d0b9b84f140b55cb75815a6aa6c983b51a4f" have entirely different histories.
c2e6420243
...
d0b4d0b9b8
@ -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
|
||||
@ -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
|
||||
@ -64,7 +64,7 @@ function plotAngularDistributionForDifferentBeta(obj, Beta)
|
||||
leg = legend;
|
||||
hXLabel = xlabel('\theta (rad)');
|
||||
hYLabel = ylabel('J(\theta)');
|
||||
hTitle = sgtitle('Angular Distribution (Transition Flow)');
|
||||
hTitle = sgtitle('Angular Distribution (Free Molecular Flow)');
|
||||
|
||||
set([hXLabel, hYLabel, leg] , ...
|
||||
'FontSize' , 14 );
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
function plotCaptureVelocityVsAngle(obj)
|
||||
|
||||
theta = linspace(0, obj.NozzleExitDivergence, 1000);
|
||||
theta = linspace(0, obj.MOTExitDivergence, 40);
|
||||
CaptureVelocity = zeros(length(theta),3);
|
||||
|
||||
for i=1:length(theta)
|
||||
@ -20,16 +20,16 @@ function plotCaptureVelocityVsAngle(obj)
|
||||
screensize = get(0,'ScreenSize');
|
||||
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)');
|
||||
hTitle = sgtitle('Capture velocity for different angles of divergence');
|
||||
hTitle = sgtitle('Capture Velocity of an atomic beam from (0,0,0)');
|
||||
|
||||
set([hXLabel, hYLabel] , ...
|
||||
'FontSize' , 14 );
|
||||
set( hTitle , ...
|
||||
'FontSize' , 14 );
|
||||
'FontSize' , 18 );
|
||||
|
||||
grid on
|
||||
Helper.bringFiguresWithTagInForeground();
|
||||
|
||||
@ -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
|
||||
|
||||
@ -17,22 +17,21 @@ function plotFreeMolecularFluxVsTemp(obj, Temperature)
|
||||
|
||||
for i=1:length(Temperature)
|
||||
beta = linspace(0.01,0.5,200);
|
||||
L = 2*obj.NozzleRadius./beta;
|
||||
obj.OvenTemperature = Temperature(i);
|
||||
flux = zeros(1,length(beta));
|
||||
for j=1:length(beta)
|
||||
obj.Beta = beta(j);
|
||||
[ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor();
|
||||
flux(j) = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux();
|
||||
flux = zeros(1,length(L));
|
||||
for j=1:length(L)
|
||||
obj.NozzleLength = L(j);
|
||||
flux(j) = obj.calculateFreeMolecularRegimeFlux();
|
||||
end
|
||||
plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5)
|
||||
end
|
||||
set(gca,'yscale','log')
|
||||
|
||||
obj.reinitializeSimulator();
|
||||
[ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor();
|
||||
|
||||
xline(obj.Beta, 'k--','Linewidth', 0.5);
|
||||
fmf = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux();
|
||||
fmf = obj.calculateFreeMolecularRegimeFlux();
|
||||
yline(fmf, 'k--', 'Linewidth', 1.5);
|
||||
textstring = [sprintf('%1.e',fmf) ' atoms/s for ' '$$ \beta $$ = ' num2str(obj.Beta, '%.2f') sprintf(' @ %.2f K', obj.OvenTemperatureinKelvin)];
|
||||
txt = text((obj.Beta + 0.05*obj.Beta), (max(fmf) + 0.2*fmf), textstring, 'Interpreter','latex');
|
||||
|
||||
@ -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
|
||||
@ -50,20 +50,20 @@ function plotPhaseSpaceWithAccelerationField(obj, MaximumVelocity, NumberOfBins,
|
||||
%-------------------------------------------------------------------------
|
||||
|
||||
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)
|
||||
x =-L/2;
|
||||
vx = Y(i)*cos(Theta);
|
||||
vz = Y(i)*sin(Theta);
|
||||
r = [x,0,z];
|
||||
v = [vx,0,vz];
|
||||
ParticleDynamicalQuantities(i,:,:) = obj.solver(r, v);
|
||||
[Trajectories(i,:,:),~] = obj.solver(r, v);
|
||||
end
|
||||
|
||||
hold on
|
||||
|
||||
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
|
||||
|
||||
hold off
|
||||
|
||||
@ -41,14 +41,14 @@ function plotPositionAndVelocitySampling(obj, NumberOfBins)
|
||||
legend('FontSize', 14)
|
||||
|
||||
subplot(3,2,4)
|
||||
histogram(initialVelocities(:, 2),NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component')
|
||||
xlabel('Velocities (m/s)','FontSize', 14)
|
||||
histogram(initialVelocities(:, 2)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component')
|
||||
xlabel('Velocities (mm/s)','FontSize', 14)
|
||||
ylabel('Counts','FontSize', 14)
|
||||
legend('FontSize', 14)
|
||||
|
||||
subplot(3,2,6)
|
||||
histogram(initialVelocities(:, 3),NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component')
|
||||
xlabel('Velocities (m/s)','FontSize', 14)
|
||||
histogram(initialVelocities(:, 3)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component')
|
||||
xlabel('Velocities (mm/s)','FontSize', 14)
|
||||
ylabel('Counts','FontSize', 14)
|
||||
legend('FontSize', 14)
|
||||
|
||||
|
||||
@ -10,10 +10,7 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
|
||||
addParameter(p, 'XLabelString', 'X parameter', @ischar)
|
||||
addParameter(p, 'ErrorsForYQuantity', false, @islogical)
|
||||
addParameter(p, 'ErrorsArray', [], @isvector)
|
||||
addParameter(p, 'CIForYQuantity', false, @islogical)
|
||||
addParameter(p, 'CIArray', [], @ismatrix)
|
||||
addParameter(p, 'RescalingFactorForYQuantity', 1, @isscalar)
|
||||
addParameter(p, 'RemoveOutliers', false, @islogical)
|
||||
addParameter(p, 'YLabelString', 'Y parameter', @ischar)
|
||||
addParameter(p, 'TitleString', 'One-Parameter Scan', @ischar)
|
||||
|
||||
@ -25,10 +22,7 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
|
||||
YQuantity = p.Results.QuantityOfInterestArray;
|
||||
ErrorsForYQuantity = p.Results.ErrorsForYQuantity;
|
||||
ErrorsArray = p.Results.ErrorsArray;
|
||||
CIForYQuantity = p.Results.CIForYQuantity;
|
||||
CIArray = p.Results.CIArray;
|
||||
RescalingFactorForYQuantity = p.Results.RescalingFactorForYQuantity;
|
||||
RemoveOutliers = p.Results.RemoveOutliers;
|
||||
YLabelString = p.Results.YLabelString;
|
||||
TitleString = p.Results.TitleString;
|
||||
|
||||
@ -45,37 +39,16 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
|
||||
screensize = get(0,'ScreenSize');
|
||||
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;
|
||||
RescaledYQuantity = YQuantity .* RescalingFactorForYQuantity;
|
||||
|
||||
hold on
|
||||
|
||||
if ErrorsForYQuantity
|
||||
RescaledErrorsArray = ErrorsArray .* RescalingFactorForYQuantity;
|
||||
errorbar(RescaledXParameter, RescaledYQuantity, RescaledErrorsArray, 'o', 'Linewidth', 1.5, 'MarkerFaceColor', '#0071BB')
|
||||
errorbar(RescaledXParameter, RescaledYQuantity, RescaledErrorsArray)
|
||||
else
|
||||
plot(RescaledXParameter, RescaledYQuantity, '--o', 'Linewidth', 1.5);
|
||||
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);
|
||||
hYLabel = ylabel(YLabelString);
|
||||
hTitle = sgtitle(TitleString);
|
||||
@ -85,6 +58,5 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
|
||||
set( hTitle , ...
|
||||
'FontSize' , 18 );
|
||||
|
||||
grid on
|
||||
Helper.bringFiguresWithTagInForeground();
|
||||
end
|
||||
@ -24,18 +24,18 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
BlueDetuning;
|
||||
BlueBeamRadius;
|
||||
BlueBeamWaist;
|
||||
BlueWaveNumber;
|
||||
BlueWaveVector;
|
||||
BlueSaturationIntensity;
|
||||
|
||||
OrangePower;
|
||||
OrangeDetuning;
|
||||
OrangeBeamRadius;
|
||||
OrangeBeamWaist;
|
||||
OrangeWaveNumber;
|
||||
OrangeWaveVector;
|
||||
OrangeSaturationIntensity;
|
||||
|
||||
CoolingBeamPower;
|
||||
CoolingBeamWaveNumber;
|
||||
CoolingBeamWaveVector;
|
||||
CoolingBeamLinewidth;
|
||||
CoolingBeamDetuning;
|
||||
CoolingBeamRadius;
|
||||
@ -49,7 +49,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
SidebandBeamSaturationIntensity;
|
||||
|
||||
PushBeamPower;
|
||||
PushBeamWaveNumber;
|
||||
PushBeamWaveVector;
|
||||
PushBeamLinewidth;
|
||||
PushBeamDetuning;
|
||||
PushBeamRadius;
|
||||
@ -71,16 +71,16 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
CaptureVelocity;
|
||||
VelocityCutoff;
|
||||
ClausingFactor;
|
||||
ReducedClausingFactor;
|
||||
ReducedFlux;
|
||||
TimeSpentInInteractionRegion;
|
||||
ThetaArray;
|
||||
AngularDistribution;
|
||||
NormalizationConstantForAngularDistribution;
|
||||
|
||||
%Flags
|
||||
SpontaneousEmission;
|
||||
Sideband;
|
||||
ZeemanSlowerBeam;
|
||||
Gravity;
|
||||
BackgroundCollision;
|
||||
AtomicBeamCollision;
|
||||
|
||||
DebugMode;
|
||||
DoSave;
|
||||
@ -124,7 +124,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
@islogical);
|
||||
addParameter(p, 'Gravity', false,...
|
||||
@islogical);
|
||||
addParameter(p, 'BackgroundCollision', false,...
|
||||
addParameter(p, 'AtomicBeamCollision', true,...
|
||||
@islogical);
|
||||
addParameter(p, 'DebugMode', false,...
|
||||
@islogical);
|
||||
@ -144,7 +144,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
s.Sideband = p.Results.Sideband;
|
||||
s.ZeemanSlowerBeam = p.Results.ZeemanSlowerBeam;
|
||||
s.Gravity = p.Results.Gravity;
|
||||
s.BackgroundCollision = p.Results.BackgroundCollision;
|
||||
s.AtomicBeamCollision = p.Results.AtomicBeamCollision;
|
||||
|
||||
s.DebugMode = p.Results.DebugMode;
|
||||
s.DoSave = p.Results.SaveData;
|
||||
@ -177,7 +177,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
ret = this.SimulationTime;
|
||||
end
|
||||
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;
|
||||
end
|
||||
function ret = get.NumberOfAtoms(this)
|
||||
@ -282,11 +282,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
function ret = get.BlueBeamWaist(this)
|
||||
ret = this.BlueBeamWaist;
|
||||
end
|
||||
function set.BlueWaveNumber(this, val)
|
||||
this.BlueWaveNumber = val;
|
||||
function set.BlueWaveVector(this, val)
|
||||
this.BlueWaveVector = val;
|
||||
end
|
||||
function ret = get.BlueWaveNumber(this)
|
||||
ret = this.BlueWaveNumber;
|
||||
function ret = get.BlueWaveVector(this)
|
||||
ret = this.BlueWaveVector;
|
||||
end
|
||||
function set.BlueSaturationIntensity(this, val)
|
||||
this.BlueSaturationIntensity = val;
|
||||
@ -319,11 +319,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
function ret = get.OrangeBeamWaist(this)
|
||||
ret = this.OrangeBeamWaist;
|
||||
end
|
||||
function set.OrangeWaveNumber(this, val)
|
||||
this.OrangeWaveNumber = val;
|
||||
function set.OrangeWaveVector(this, val)
|
||||
this.OrangeWaveVector = val;
|
||||
end
|
||||
function ret = get.OrangeWaveNumber(this)
|
||||
ret = this.OrangeWaveNumber;
|
||||
function ret = get.OrangeWaveVector(this)
|
||||
ret = this.OrangeWaveVector;
|
||||
end
|
||||
function set.OrangeSaturationIntensity(this, val)
|
||||
this.OrangeSaturationIntensity = val;
|
||||
@ -356,11 +356,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
function ret = get.CoolingBeamWaist(this)
|
||||
ret = this.CoolingBeamWaist;
|
||||
end
|
||||
function set.CoolingBeamWaveNumber(this, val)
|
||||
this.CoolingBeamWaveNumber = val;
|
||||
function set.CoolingBeamWaveVector(this, val)
|
||||
this.CoolingBeamWaveVector = val;
|
||||
end
|
||||
function ret = get.CoolingBeamWaveNumber(this)
|
||||
ret = this.CoolingBeamWaveNumber;
|
||||
function ret = get.CoolingBeamWaveVector(this)
|
||||
ret = this.CoolingBeamWaveVector;
|
||||
end
|
||||
function set.CoolingBeamLinewidth(this, val)
|
||||
this.CoolingBeamLinewidth = val;
|
||||
@ -430,11 +430,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
function ret = get.PushBeamWaist(this)
|
||||
ret = this.PushBeamWaist;
|
||||
end
|
||||
function set.PushBeamWaveNumber(this, val)
|
||||
this.PushBeamWaveNumber= val;
|
||||
function set.PushBeamWaveVector(this, val)
|
||||
this.PushBeamWaveVector = val;
|
||||
end
|
||||
function ret = get.PushBeamWaveNumber(this)
|
||||
ret = this.PushBeamWaveNumber;
|
||||
function ret = get.PushBeamWaveVector(this)
|
||||
ret = this.PushBeamWaveVector;
|
||||
end
|
||||
function set.PushBeamLinewidth(this, val)
|
||||
this.PushBeamLinewidth = val;
|
||||
@ -529,23 +529,31 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
function ret = get.ClausingFactor(this)
|
||||
ret = this.ClausingFactor;
|
||||
end
|
||||
function set.ReducedClausingFactor(this,val)
|
||||
this.ReducedClausingFactor = val;
|
||||
function set.AngularDistribution(this,val)
|
||||
this.AngularDistribution = val;
|
||||
end
|
||||
function ret = get.ReducedClausingFactor(this)
|
||||
ret = this.ReducedClausingFactor;
|
||||
function ret = get.AngularDistribution(this)
|
||||
ret = this.AngularDistribution;
|
||||
end
|
||||
function set.ReducedFlux(this,val)
|
||||
this.ReducedFlux = val;
|
||||
function set.ThetaArray(this,val)
|
||||
this.ThetaArray = val;
|
||||
end
|
||||
function ret = get.ReducedFlux(this)
|
||||
ret = this.ReducedFlux;
|
||||
function ret = get.ThetaArray(this)
|
||||
ret = this.ThetaArray;
|
||||
end
|
||||
function set.TimeSpentInInteractionRegion(this,val)
|
||||
this.TimeSpentInInteractionRegion = val;
|
||||
function set.NormalizationConstantForAngularDistribution(this,val)
|
||||
this.NormalizationConstantForAngularDistribution = val;
|
||||
end
|
||||
function ret = get.TimeSpentInInteractionRegion(this)
|
||||
ret = this.TimeSpentInInteractionRegion;
|
||||
function ret = get.NormalizationConstantForAngularDistribution(this)
|
||||
ret = this.NormalizationConstantForAngularDistribution;
|
||||
end
|
||||
|
||||
function set.AtomicBeamCollision(this,val)
|
||||
this.AtomicBeamCollision=val;
|
||||
end
|
||||
|
||||
function ret=get.AtomicBeamCollision(this)
|
||||
ret=this.AtomicBeamCollision;
|
||||
end
|
||||
|
||||
function set.DebugMode(this, val)
|
||||
@ -578,19 +586,19 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
|
||||
methods
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
function ret = get.OvenTemperatureinKelvin(this)
|
||||
@ -599,7 +607,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
|
||||
|
||||
function ret = get.AverageVelocity(this)
|
||||
%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
|
||||
|
||||
function ret = get.AtomicBeamDensity(this)
|
||||
|
||||
@ -4,17 +4,17 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
|
||||
WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint);
|
||||
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;
|
||||
|
||||
s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2));
|
||||
a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* (s_push/(1+s_push));
|
||||
a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveVector * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* ...
|
||||
(SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2)));
|
||||
|
||||
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
|
||||
a_scatter = [0,0,0];
|
||||
end
|
||||
|
||||
@ -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 = 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);
|
||||
|
||||
if NumberOfScatteringEvents > 0
|
||||
ret = Vector.*((Helper.PhysicsConstants.PlanckConstantReduced * WaveNumber) / ...
|
||||
ret = Vector.*((Helper.PhysicsConstants.PlanckConstantReduced * WaveVector) / ...
|
||||
(Helper.PhysicsConstants.Dy164Mass * this.TimeStep)).* sqrt(NumberOfScatteringEvents);
|
||||
else
|
||||
ret = zeros(1,3);
|
||||
|
||||
@ -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
|
||||
@ -1,15 +1,14 @@
|
||||
function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
|
||||
|
||||
switch this.SimulationMode
|
||||
case "2D"
|
||||
VelocityUnitVector = VelocityVector./norm(VelocityVector);
|
||||
VelocityVector = VelocityVector./norm(VelocityVector);
|
||||
UpperLimit = 500;
|
||||
LowerLimit = 0;
|
||||
this.AtomicBeamCollision = false;
|
||||
|
||||
for Index = 1:500
|
||||
InitialVelocity = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector;
|
||||
ParticleDynamicalQuantities = this.solver(PositionVector, InitialVelocity);
|
||||
FinalPositionVector = ParticleDynamicalQuantities(end, 1:3);
|
||||
InitialVelocity = 0.5 * (UpperLimit + LowerLimit) * VelocityVector;
|
||||
[~, FinalDynamicalQuantities] = this.solver(PositionVector, InitialVelocity);
|
||||
FinalPositionVector = FinalDynamicalQuantities(1:3);
|
||||
if rssq(FinalPositionVector) <= this.OvenDistance
|
||||
LowerLimit = 0.5 * (UpperLimit + LowerLimit);
|
||||
else
|
||||
@ -17,12 +16,9 @@ function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
|
||||
end
|
||||
|
||||
if UpperLimit - LowerLimit < 1
|
||||
ret = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector;
|
||||
ret = InitialVelocity;
|
||||
break;
|
||||
end
|
||||
end
|
||||
case "3D"
|
||||
% Development In progress
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
@ -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
|
||||
@ -3,8 +3,6 @@ function ret = calculateFreeMolecularRegimeFlux(this)
|
||||
%See Expected atomic flux section in Barbiero
|
||||
Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature
|
||||
Dy164DensityinOven = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin);
|
||||
ret = Helper.PhysicsConstants.Dy164IsotopicAbundance * 1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2;
|
||||
% Needs to be multiplied with the "Clausing Factor" which here would be
|
||||
% the probability not for the full solid angle but the angle subtended
|
||||
% by the aperture of the oven at its mouth.
|
||||
ClausingFactorApproximation = (8 * this.NozzleRadius) / (3 * this.NozzleLength);
|
||||
ret = Helper.PhysicsConstants.Dy164IsotopicAbundance * 1/4 * Dy164DensityinOven * this.AverageVelocity * ClausingFactorApproximation * pi * this.NozzleRadius.^2;
|
||||
end
|
||||
@ -1,38 +1,84 @@
|
||||
function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities)
|
||||
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);
|
||||
function [LoadingRate, StandardError] = calculateLoadingRate(this, FinalDynamicalQuantities)
|
||||
|
||||
% Include the stochastic process of background collisions
|
||||
for AtomIndex = 1:n
|
||||
TimeCounts(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
|
||||
end
|
||||
this.TimeSpentInInteractionRegion = mean(TimeCounts);
|
||||
for AtomIndex = 1:n
|
||||
CollisionEvents(AtomIndex) = this.computeCollisionProbability();
|
||||
NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms);
|
||||
AutocorrelationFunction = zeros(1, this.NumberOfAtoms);
|
||||
|
||||
for i = 1:this.NumberOfAtoms
|
||||
FinalPosition = FinalDynamicalQuantities(i,1:3);
|
||||
DivergenceAngle = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2)));
|
||||
if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0)
|
||||
if i == 1
|
||||
NumberOfLoadedAtoms(1) = 1;
|
||||
else
|
||||
NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1) + 1;
|
||||
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(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
|
||||
Velocity = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
|
||||
if this.exitCondition(Position, Velocity, CollisionEvents(AtomIndex))
|
||||
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
|
||||
else
|
||||
if i == 1
|
||||
NumberOfLoadedAtoms(1) = 0;
|
||||
else
|
||||
NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1);
|
||||
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"
|
||||
% Development In progress
|
||||
if AutocorrelationFunction(1)~=0 % To handle cases where there is no atom loading
|
||||
|
||||
AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1);
|
||||
x = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms);
|
||||
[FitObject,~] = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100);
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
@ -1,6 +1,12 @@
|
||||
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
|
||||
ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2);
|
||||
|
||||
@ -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,9 +28,9 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
|
||||
|
||||
B = sign(LocalMagneticField(1:3) * WaveVectorEndPoint(i,1:3)') * LocalMagneticField(4);
|
||||
|
||||
ZeemanShift = this.LandegFactor * this.MagneticSubLevel * (Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced) * B;
|
||||
ZeemanShift = this.LandegFactor * this.MagneticSubLevel * Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced * B;
|
||||
|
||||
DopplerShift = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveNumber;
|
||||
DopplerShift = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveVector;
|
||||
|
||||
Delta_Cooling(i*2-1) = 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];
|
||||
|
||||
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) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Cooling(2*i) /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)^2 / this.CoolingBeamLinewidth^2);
|
||||
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+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2)) /(1 + 4 * (Delta_Sideband(2*i)/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)^2 / this.CoolingBeamLinewidth^2);
|
||||
end
|
||||
end
|
||||
|
||||
@ -56,28 +56,28 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
|
||||
|
||||
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));
|
||||
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));
|
||||
|
||||
if this.SpontaneousEmission
|
||||
a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ...
|
||||
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, Delta_Cooling(2*i), 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.CoolingBeamWaveVector);
|
||||
else
|
||||
a_scattering = [0,0,0];
|
||||
end
|
||||
|
||||
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));
|
||||
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));
|
||||
|
||||
if this.SpontaneousEmission
|
||||
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+4), TotalSaturationParameter, Delta_Cooling(2*i), 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.CoolingBeamWaveVector);
|
||||
else
|
||||
a_scattering = [0,0,0];
|
||||
end
|
||||
|
||||
@ -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
|
||||
@ -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
|
||||
@ -1,4 +1,4 @@
|
||||
function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOneParameterScan(this, ParameterName, ParameterArray, varargin)
|
||||
function [LoadingRateArray, StandardErrorArray] = doOneParameterScan(this, ParameterName, ParameterArray, varargin)
|
||||
|
||||
p = inputParser;
|
||||
p.KeepUnmatched = true;
|
||||
@ -30,7 +30,6 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOne
|
||||
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)));
|
||||
@ -45,18 +44,15 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOne
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
this.setInitialConditions(options{:});
|
||||
tic
|
||||
[LoadingRate, StandardError, ConfidenceInterval] = this.runSimulation();
|
||||
[LoadingRate, StandardError] = 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')];
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ...
|
||||
function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ...
|
||||
SecondParameterName, SecondParameterArray, varargin)
|
||||
|
||||
p = inputParser;
|
||||
@ -38,7 +38,6 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwo
|
||||
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)));
|
||||
@ -59,11 +58,9 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwo
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
this.setInitialConditions(options{:});
|
||||
tic
|
||||
[LoadingRate, StandardError, ConfidenceInterval] = this.runSimulation();
|
||||
[LoadingRate, StandardError] = this.runSimulation();
|
||||
LoadingRateArray(i,j) = LoadingRate;
|
||||
StandardErrorArray(i,j) = StandardError;
|
||||
ConfidenceIntervalArray(i, j, 1) = ConfidenceInterval(1);
|
||||
ConfidenceIntervalArray(i, j, 2) = ConfidenceInterval(2);
|
||||
end
|
||||
end
|
||||
|
||||
@ -71,7 +68,6 @@ function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwo
|
||||
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')];
|
||||
|
||||
@ -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
|
||||
@ -1,12 +1,7 @@
|
||||
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
|
||||
|
||||
|
||||
@ -1,25 +1,6 @@
|
||||
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)/integral(VelocityDistribution,0,Inf);
|
||||
|
||||
this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
|
||||
|
||||
ret = zeros(n,3);
|
||||
SampledVelocityMagnitude = zeros(n,1);
|
||||
SampledPolarAngle = zeros(n,1);
|
||||
@ -32,25 +13,27 @@ function ret = initialVelocitySampling(this)
|
||||
else
|
||||
MaximumVelocityAllowed = MostProbableVelocity;
|
||||
end
|
||||
NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed);
|
||||
ProbabilityOfMaximumVelocityAllowed = this.velocityDistributionFunction(MaximumVelocityAllowed);
|
||||
ProbabilityOfMaximumDivergenceAngleAllowed = 0.98 * this.NormalizationConstantForAngularDistribution * max(this.AngularDistribution .* sin(this.ThetaArray));
|
||||
|
||||
parfor i = 1:n
|
||||
% Rejection Sampling of speed
|
||||
y = rand(1);
|
||||
y = ProbabilityOfMaximumVelocityAllowed * 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);
|
||||
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);
|
||||
z = this.MOTExitDivergence * 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
|
||||
w = rand(1);
|
||||
z = this.NozzleExitDivergence * rand(1);
|
||||
while w > (this.NormalizationConstantForAngularDistribution * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value
|
||||
z = this.MOTExitDivergence * rand(1);
|
||||
w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
|
||||
end
|
||||
SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
|
||||
|
||||
@ -60,8 +43,5 @@ function ret = initialVelocitySampling(this)
|
||||
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
|
||||
|
||||
|
||||
@ -3,9 +3,8 @@ function reinitializeSimulator(this)
|
||||
pc = Helper.PhysicsConstants;
|
||||
%% SIMULATION PARAMETERS
|
||||
this.NozzleLength = 60e-3;
|
||||
this.NozzleRadius = 2.60e-3;
|
||||
this.NozzleRadius = 2.50e-3;
|
||||
this.Beta = 2 * (this.NozzleRadius/this.NozzleLength);
|
||||
this.ClausingFactor = this.calculateClausingFactor();
|
||||
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
|
||||
@ -13,25 +12,20 @@ function reinitializeSimulator(this)
|
||||
% 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.MOTDistance = 320e-3; % Distance between the 2-D MOT the 3-D MOT
|
||||
this.BlueWaveVector = 2*pi/pc.BlueWavelength;
|
||||
this.BlueSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.BlueLinewidth/3/(pc.BlueWavelength)^3/10;
|
||||
this.OrangeWaveVector = 2*pi/pc.OrangeWavelength;
|
||||
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
|
||||
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.TimeSpentInInteractionRegion = this.SimulationTime;
|
||||
this.TotalPower = 0.8;
|
||||
this.OrangeBeamRadius = 1.2e-3;
|
||||
this.PushBeamRadius = 1.2e-3;
|
||||
this.MOTExitDivergence = 0.016; % The limitation angle between 2D-MOT and 3D-MOT
|
||||
this.TotalPower = 0.4;
|
||||
this.OrangeBeamRadius = 1.2e-03;
|
||||
this.PushBeamRadius = 0.000;
|
||||
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.ZeemanSlowerBeamRadius = 1;
|
||||
end
|
||||
@ -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
|
||||
this.InitialPositions = this.initialPositionSampling();
|
||||
% - sampling the velocity distribution
|
||||
@ -8,18 +38,18 @@ function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this)
|
||||
|
||||
%% Solve ODE
|
||||
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
|
||||
ParticleDynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
|
||||
FinalDynamicalQuantities = zeros(n,9);
|
||||
Positions = this.InitialPositions;
|
||||
Velocities = this.InitialVelocities;
|
||||
parfor Index = 1:this.NumberOfAtoms
|
||||
ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
|
||||
parfor Index = 1:n
|
||||
ret = this.solver(Positions(Index,:), Velocities(Index,:));
|
||||
FinalDynamicalQuantities(Index,:) = ret(end,:);
|
||||
progressbar.PB_iterate();
|
||||
end
|
||||
clear Index
|
||||
|
||||
%% Calculate the Loading Rate
|
||||
[LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ParticleDynamicalQuantities);
|
||||
[LoadingRate, StandardError] = this.calculateLoadingRate(FinalDynamicalQuantities);
|
||||
end
|
||||
@ -4,27 +4,27 @@ function setInitialConditions(this, varargin)
|
||||
p.KeepUnmatched = true;
|
||||
addParameter(p, 'NumberOfAtoms', 5000,...
|
||||
@(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)));
|
||||
addParameter(p, 'BlueDetuning', -1.39*Helper.PhysicsConstants.BlueLinewidth,...
|
||||
addParameter(p, 'BlueDetuning', -1.5*Helper.PhysicsConstants.BlueLinewidth,...
|
||||
@(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)));
|
||||
addParameter(p, 'SidebandPower', 0.5*this.TotalPower,...
|
||||
addParameter(p, 'SidebandPower', 200e-3,...
|
||||
@(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)));
|
||||
addParameter(p, 'SidebandBeamWaist', 15e-3,...
|
||||
addParameter(p, 'SidebandBeamWaist', 12e-3,...
|
||||
@(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)));
|
||||
addParameter(p, 'PushBeamDetuning', -5*Helper.PhysicsConstants.OrangeLinewidth,...
|
||||
addParameter(p, 'PushBeamDetuning', 0,...
|
||||
@(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,...
|
||||
addParameter(p, 'OrangeDetuning', 12e-3,...
|
||||
@(x) assert(isnumeric(x) && isscalar(x)));
|
||||
addParameter(p, 'OrangeBeamWaist', 12e-3,...
|
||||
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
||||
@ -56,19 +56,24 @@ function setInitialConditions(this, varargin)
|
||||
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.CoolingBeamWaveVector = this.BlueWaveVector;
|
||||
this.CoolingBeamDetuning = this.BlueDetuning;
|
||||
this.CoolingBeamRadius = this.BlueBeamRadius;
|
||||
this.CoolingBeamWaist = this.BlueBeamWaist;
|
||||
this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity;
|
||||
this.SidebandBeamRadius = this.BlueBeamRadius;
|
||||
this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity;
|
||||
this.PushBeamLinewidth = Helper.PhysicsConstants.OrangeLinewidth;
|
||||
this.PushBeamWaveVector = this.OrangeWaveVector;
|
||||
this.PushBeamDetuning = this.OrangeDetuning;
|
||||
this.PushBeamSaturationIntensity = this.OrangeSaturationIntensity;
|
||||
this.LandegFactor = Helper.PhysicsConstants.BlueLandegFactor;
|
||||
this.MagneticSubLevel = 1;
|
||||
case "3D"
|
||||
|
||||
@ -1,16 +1,22 @@
|
||||
function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelocity)
|
||||
|
||||
function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity)
|
||||
if this.Gravity
|
||||
g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration];
|
||||
else
|
||||
g = 0;
|
||||
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)
|
||||
|
||||
ParticleDynamicalQuantities(i,1:3) = InitialPosition;
|
||||
ParticleDynamicalQuantities(i,4:6) = InitialVelocity;
|
||||
ParticleTrajectory(i,1:3) = InitialPosition;
|
||||
ParticleTrajectory(i,4:6) = InitialVelocity;
|
||||
|
||||
rt = InitialPosition;
|
||||
vt = InitialVelocity;
|
||||
@ -35,7 +41,14 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelo
|
||||
|
||||
InitialPosition = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./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
|
||||
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
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 ...
|
||||
* this.OvenTemperatureinKelvin)));
|
||||
end
|
||||
@ -7,10 +7,12 @@ OptionsStruct.SimulationMode = '2D';
|
||||
OptionsStruct.TimeStep = 50e-06; % in s
|
||||
OptionsStruct.SimulationTime = 4e-03; % in s
|
||||
OptionsStruct.SpontaneousEmission = true;
|
||||
OptionsStruct.Sideband = false;
|
||||
OptionsStruct.Sideband = true;
|
||||
OptionsStruct.ZeemanSlowerBeam = false;
|
||||
OptionsStruct.Gravity = true;
|
||||
OptionsStruct.BackgroundCollision = true;
|
||||
OptionsStruct.SaveData = false;
|
||||
OptionsStruct.AtomicBeamCollision = true;
|
||||
OptionsStruct.DebugMode = false;
|
||||
OptionsStruct.SaveData = true;
|
||||
OptionsStruct.SaveDirectory = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation';
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
|
||||
@ -40,18 +42,9 @@ clear OptionsStruct
|
||||
[LoadingRate, ~] = Simulator.runSimulation();
|
||||
|
||||
%% - Plot initial distribution
|
||||
Simulator.setInitialConditions();
|
||||
% - sampling the position distribution
|
||||
Simulator.InitialPositions = Simulator.initialPositionSampling();
|
||||
% - sampling the velocity distribution
|
||||
Simulator.InitialVelocities = Simulator.initialVelocitySampling();
|
||||
NumberOfBins = 100;
|
||||
Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins);
|
||||
|
||||
%% - Plot distributions of magnitude and direction of initial velocities
|
||||
NumberOfBins = 50;
|
||||
Plotting.plotInitialVeloctiySamplingVsAngle(Simulator, NumberOfBins)
|
||||
|
||||
%% - Plot Magnetic Field
|
||||
XAxisRange = [-5 5];
|
||||
YAxisRange = [-5 5];
|
||||
@ -63,7 +56,7 @@ TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius
|
||||
Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
|
||||
|
||||
%% - Plot the Free Molecular Flux for different temperatures
|
||||
Temperature = [950, 1000, 1050]; % Temperature
|
||||
Temperature = [900, 950, 1000, 1050, 1100]; % Temperature'
|
||||
Plotting.plotFreeMolecularFluxVsTemp(Simulator,Temperature)
|
||||
|
||||
%% - 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)
|
||||
|
||||
%% - Plot Capture Velocity
|
||||
|
||||
Simulator.setInitialConditions();
|
||||
Simulator.SimulationTime = 15*10^(-4);
|
||||
Simulator.TimeStep = 1.5*10^(-6);
|
||||
Simulator.MOTExitDivergence = 15/360*2*pi;
|
||||
Simulator.MagneticGradient = 0.4;
|
||||
|
||||
Plotting.plotCaptureVelocityVsAngle(Simulator);
|
||||
|
||||
%% - Plot Phase Space with Acceleration Field
|
||||
|
||||
Simulator.NumberOfAtoms = 200;
|
||||
Simulator.NumberOfAtoms = 100;
|
||||
MaximumVelocity = 150;
|
||||
NumberOfBins = 200; %Along each axis
|
||||
IncidentAtomDirection = 0*2*pi/360;
|
||||
@ -87,7 +86,7 @@ Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberO
|
||||
|
||||
% ONE-PARAMETER SCAN
|
||||
|
||||
NumberOfPointsForParam = 10; %iterations of the simulation
|
||||
NumberOfPointsForParam = 20; %iterations of the simulation
|
||||
% Scan Cooling Beam Power
|
||||
PowerArray = linspace(0.1, 1.0, NumberOfPointsForParam) * Simulator.TotalPower;
|
||||
% Scan Cooling Beam Detuning
|
||||
@ -96,12 +95,11 @@ PowerArray = linspace(0.1, 1.0, NumberOfPointsForParam) * Simulator.TotalPow
|
||||
OptionsStruct = struct;
|
||||
OptionsStruct.ChangeInitialConditions = true;
|
||||
OptionsStruct.ParameterNameArray = {'NumberOfAtoms'};
|
||||
OptionsStruct.ParameterValueArray = {5000};
|
||||
OptionsStruct.ParameterValueArray = {10000};
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
|
||||
|
||||
tStart = tic;
|
||||
[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
|
||||
[LoadingRateArray, ~] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
|
||||
tEnd = toc(tStart);
|
||||
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
|
||||
|
||||
@ -109,19 +107,15 @@ clear OptionsStruct
|
||||
|
||||
% - Plot results
|
||||
|
||||
ParameterArray = PowerArray;
|
||||
ParameterArray = PowerArray .* Simulator.TotalPower;
|
||||
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.RescalingFactorForYQuantity = 1e-16;
|
||||
OptionsStruct.ErrorsForYQuantity = false;
|
||||
OptionsStruct.YLabelString = 'Loading rate (x 10^{16} atoms/s)';
|
||||
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
|
||||
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
@ -132,11 +126,11 @@ clear OptionsStruct
|
||||
|
||||
%% TWO-PARAMETER SCAN
|
||||
|
||||
NumberOfPointsForParam = 10; %iterations of the simulation
|
||||
NumberOfPointsForParam = 20; %iterations of the simulation
|
||||
NumberOfPointsForSecondParam = 10;
|
||||
|
||||
% 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;
|
||||
BluePowerArray = Simulator.TotalPower - SidebandPowerArray;
|
||||
|
||||
@ -147,11 +141,11 @@ OptionsStruct.RelatedParameterName = 'BluePower';
|
||||
OptionsStruct.RelatedParameterArray = BluePowerArray;
|
||||
OptionsStruct.ChangeInitialConditions = true;
|
||||
OptionsStruct.ParameterNameArray = {'NumberOfAtoms'};
|
||||
OptionsStruct.ParameterValueArray = {5000};
|
||||
OptionsStruct.ParameterValueArray = {10000};
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
|
||||
tStart = tic;
|
||||
[LoadingRateArray, StandardErrorArray, ConfidenceInterval] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
|
||||
[LoadingRateArray, StandardErrorArray] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
|
||||
tEnd = toc(tStart);
|
||||
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
|
||||
|
||||
@ -168,8 +162,8 @@ OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.B
|
||||
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.RescalingFactorForQuantityOfInterest = 1e-16;
|
||||
OptionsStruct.ZLabelString = 'Loading rate (x 10^{16} atoms/s)';
|
||||
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
|
||||
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user