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

Compare commits

base: karthik:149d59f7c446a9f90d2a6f438297b13602b84127
Branches Tags
karthik:master
..
compare: karthik:de683c7d6b23129e7efd4ea90667d2d2f22a19fd
Branches Tags
karthik:master

No commits in common. "149d59f7c446a9f90d2a6f438297b13602b84127" and "de683c7d6b23129e7efd4ea90667d2d2f22a19fd" have entirely different histories.
149d59f7c4 ... de683c7d6b

16 changed files with 61 additions and 180 deletions
Show Stats Expand all files Collapse all files

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

@ -2,7 +2,7 @@ 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 = norm(CrossProduct) / norm(p12);
ret = rssq(CrossProduct) / rssq(p12);
%Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
%Area = One side of parallelogram X Base

107
MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m
View File

@ -1,107 +0,0 @@
function plotDynamicalQuantities(obj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition)
f_h = Helper.getFigureByTag('Trajectories Plot');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Trajectories Plot';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/50 screensize(4)/10] 1870 850];
N = obj.NumberOfAtoms;
Theta = IncidentAtomDirection;
z = IncidentAtomPosition;
obj.setInitialConditions();
L = obj.OvenDistance * 2;
T = obj.SimulationTime;
tau = obj.TimeStep;
Y = linspace(0,MaximumVelocity,N);
DynamicalQuantities = zeros(length(Y),int64(T/tau),6);
for i=1:length(Y)
x =-L/2;
vx = Y(i)*cos(Theta);
vz = Y(i)*sin(Theta);
r = [x,0,z];
v = [vx,0,vz];
DynamicalQuantities(i,:,:) = obj.solver(r, v);
end
t = linspace(0, T, T/tau) * 1e+3;
for i=1:length(Y)
subplot(3, 3, 1)
hold on
plot(t, DynamicalQuantities(i,:,1) * 1e+3,'r','linewidth',1.3)
hXLabel_1 = xlabel('Time (ms)');
hYLabel_1 = ylabel('x (mm)');
hold off
subplot(3, 3, 2)
hold on
plot(t, DynamicalQuantities(i,:,2) * 1e+3,'g','linewidth',1.3)
hXLabel_2 = xlabel('Time (ms)');
hYLabel_2 = ylabel('y (mm)');
hold off
subplot(3, 3, 3)
hold on
plot(t, DynamicalQuantities(i,:,3) * 1e+3,'b','linewidth',1.3)
hXLabel_3 = xlabel('Time (ms)');
hYLabel_3 = ylabel('z (m/s)');
hold off
subplot(3, 3, 4)
hold on
plot(t, DynamicalQuantities(i,:,4),'r','linewidth',1.3)
hXLabel_4 = xlabel('Time (ms)');
hYLabel_4 = ylabel('v_x (m/s)');
hold off
subplot(3, 3, 5)
hold on
plot(t, DynamicalQuantities(i,:,5),'g','linewidth',1.3)
hXLabel_5 = xlabel('Time (ms)');
hYLabel_5 = ylabel('v_y (m/s)');
hold off
subplot(3, 3, 6)
hold on
plot(t, DynamicalQuantities(i,:,6),'b','linewidth',1.3)
hXLabel_6 = xlabel('Time (ms)');
hYLabel_6 = ylabel('v_z (m/s)');
hold off
subplot(3, 3, 7)
hold on
plot(DynamicalQuantities(i,:,1), DynamicalQuantities(i,:,4),'r','linewidth',1.3)
hXLabel_7 = xlabel('x (mm)');
hYLabel_7 = ylabel('v_x (m/s)');
xlim([-L/2 L/2])
hold off
subplot(3, 3, 8)
hold on
plot(DynamicalQuantities(i,:,2), DynamicalQuantities(i,:,5),'g','linewidth',1.3)
hXLabel_8 = xlabel('y (mm)');
hYLabel_8 = ylabel('v_y (m/s)');
xlim([-1 1] * 1e-04)
hold off
subplot(3, 3, 9)
hold on
plot(DynamicalQuantities(i,:,3), DynamicalQuantities(i,:,6),'b','linewidth',1.3)
hXLabel_9 = xlabel('z (mm)');
hYLabel_9 = ylabel('v_z (m/s)');
xlim([-1 1] * 1e-04)
hold off
end
hold off
hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient));
set([hXLabel_1, hXLabel_2, hXLabel_3, hXLabel_4, hXLabel_5, hXLabel_6, hXLabel_7, hXLabel_8, hXLabel_9,...
hYLabel_1, hYLabel_2, hYLabel_3, hYLabel_4, hYLabel_5, hYLabel_6, hYLabel_7, hYLabel_8, hYLabel_9], ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

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

@ -49,7 +49,7 @@ function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins)
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('|v| (m/s)');
hXLabel_1 = xlabel('Magnitudes (m/s)');
hYLabel_1 = ylabel('Counts');
hold off
@ -59,7 +59,7 @@ function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins)
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('\theta (mrad)');
hXLabel_2 = xlabel(' Directions (mrad)');
hYLabel_2 = ylabel('Counts');
hold off

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

@ -1,4 +1,4 @@
function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
function plotPhaseSpaceWithAccelerationField(obj, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
f_h = Helper.getFigureByTag('Phase Space Plot');
set(groot,'CurrentFigure',f_h);
@ -49,21 +49,21 @@ function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVeloci
%-------------------------------------------------------------------------
Y = linspace(MinimumVelocity, MaximumVelocity,N);
DynamicalQuantities = zeros(length(Y),int64(T/tau),6);
Y = linspace(0,MaximumVelocity,N);
ParticleDynamicalQuantities = zeros(length(Y),int64(T/tau),6);
for i=1:length(Y)
x =-L/2;
vx = Y(i)*cos(Theta);
vz = Y(i)*sin(Theta);
r = [x,0,z];
v = [vx,0,vz];
DynamicalQuantities(i,:,:) = obj.solver(r, v);
ParticleDynamicalQuantities(i,:,:) = obj.solver(r, v);
end
hold on
for i=1:length(Y)
plot(DynamicalQuantities(i,:,1),DynamicalQuantities(i,:,4),'w','linewidth',1.3)
plot(ParticleDynamicalQuantities(i,:,1),ParticleDynamicalQuantities(i,:,4),'w','linewidth',1.3)
end
hold off

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

@ -1,4 +1,4 @@
function plotPositionAndVelocitySampling(NumberOfBins, initialPositions, initialVelocities)
function plotPositionAndVelocitySampling(obj, NumberOfBins)
f_h = Helper.getFigureByTag('RejectionSampling');
set(groot,'CurrentFigure',f_h);
@ -13,6 +13,9 @@ function plotPositionAndVelocitySampling(NumberOfBins, initialPositions, initial
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/7 screensize(4)/7] 1.357e+03 770];
initialPositions = obj.InitialPositions;
initialVelocities = obj.InitialVelocities;
subplot(3,2,1)
histogram(initialPositions(:, 1)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component')
xlabel('Positions (mm)','FontSize', 14)

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

@ -6,7 +6,8 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
SimulationTime;
NumberOfAtoms;
ParticleDynamicalQuantities;
InitialPositions;
InitialVelocities;
NozzleLength;
NozzleRadius;
@ -187,11 +188,17 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
ret = this.NumberOfAtoms;
end
function set.ParticleDynamicalQuantities(this,val)
this.ParticleDynamicalQuantities = val;
function set.InitialPositions(this,val)
this.InitialPositions = val;
end
function ret = get.ParticleDynamicalQuantities(this)
ret = this.ParticleDynamicalQuantities;
function ret = get.InitialPositions(this)
ret = this.InitialPositions;
end
function set.InitialVelocities(this,val)
this.InitialVelocities = val;
end
function ret = get.InitialVelocities(this)
ret = this.InitialVelocities;
end
function set.NozzleLength(this,val)
@ -596,7 +603,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
function ret = get.AverageVelocity(this)
%See Background collision probability section in Barbiero
ret = sqrt((8 * pi *Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass));
ret = sqrt((8*Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (pi*Helper.PhysicsConstants.Dy164Mass));
end
function ret = get.AtomicBeamDensity(this)
@ -617,8 +624,6 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
end % - getters for dependent properties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods(Access = protected)
function cp = copyElement(this)
% Shallow copy object
@ -637,18 +642,4 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
end
end
methods (Static)
% Creates an Instance of Class, ensures singleton behaviour (that there
% can only be one Instance of this class
function singleObj = getInstance(varargin)
% Creates an Instance of Class, ensures singleton behaviour
persistent localObj;
if isempty(localObj) || ~isvalid(localObj)
localObj = MOTSimulator(varargin{:});
end
singleObj = localObj;
end
end
end

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

@ -6,13 +6,12 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist);
DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * this.PushBeamWaveNumber;
DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveNumber;
Detuning = this.PushBeamDetuning - DopplerShift;
s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2));
a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5);
a_push = a_sat .* (s_push/(1+s_push));
a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* (s_push/(1+s_push));
if this.SpontaneousEmission
a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, this.PushBeamLinewidth, this.PushBeamWaveNumber);

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

@ -2,7 +2,7 @@ function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationParam
Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1];
Vector = Vector./norm(Vector);
ScatteringRate = (0.5 * Linewidth) * (SaturationParameter / (1 + TotalSaturationParameter));
ScatteringRate = (0.5 * SaturationParameter * Linewidth) / (1 + TotalSaturationParameter);
NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep);
if NumberOfScatteringEvents > 0

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

@ -3,7 +3,7 @@ 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 * pi * this.NozzleRadius.^2;
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.

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

@ -1,8 +1,7 @@
function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this)
function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities)
switch this.SimulationMode
case "2D"
n = this.NumberOfAtoms;
DynamicalQuantities = this.ParticleDynamicalQuantities;
NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep);
NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps);
TimeCounts = zeros(1, n);
@ -10,7 +9,7 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
% Include the stochastic process of background collisions
for AtomIndex = 1:n
TimeCounts(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(DynamicalQuantities(AtomIndex,:,1:3)));
TimeCounts(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
end
this.TimeSpentInInteractionRegion = mean(TimeCounts);
for AtomIndex = 1:n
@ -23,8 +22,8 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
end
for AtomIndex = 1:n
Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
Velocity = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
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

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

@ -43,11 +43,11 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
SaturationParameter = [0,0,0,0,0,0,0,0];
for i = 1:2
SaturationParameter(2*i-1) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
SaturationParameter(2*i) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i) /this.CoolingBeamLinewidth)^2);
SaturationParameter(2*i-1) = CoolingBeamLocalSaturationIntensity(1) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
SaturationParameter(2*i) = CoolingBeamLocalSaturationIntensity(2) /(1 + 4 * (Delta_Cooling(2*i) /this.CoolingBeamLinewidth)^2);
if this.Sideband
SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
SaturationParameter(2*i+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(1) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
SaturationParameter(2*i+4) = SidebandLocalSaturationIntensity(2) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
end
end

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

@ -16,7 +16,7 @@ function ret = initialVelocitySampling(this)
* velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* this.OvenTemperatureinKelvin)));
c = integral(VelocityDistribution, 0, this.VelocityCutoff);
c = integral(VelocityDistribution, 0, this.VelocityCutoff)/integral(VelocityDistribution,0,Inf);
this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();

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

@ -5,6 +5,7 @@ function reinitializeSimulator(this)
this.NozzleLength = 60e-3;
this.NozzleRadius = 2.60e-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
@ -22,7 +23,8 @@ function reinitializeSimulator(this)
Theta_Aperture = 15/360*2*pi; % The limitation angle of the second aperture in the oven
this.NozzleExitDivergence = min(Theta_Nozzle,Theta_Aperture);
this.MOTExitDivergence = 16e-3; % The limitation angle between 2D-MOT and 3D-MOT
this.TotalPower = 0.6;
this.TimeSpentInInteractionRegion = this.SimulationTime;
this.TotalPower = 0.8;
this.OrangeBeamRadius = 1.2e-3;
this.PushBeamRadius = 1.2e-3;
this.PushBeamDistance = 0.32;

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

@ -2,23 +2,24 @@ function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this)
%% - Sampling for initial positions and velocities
% - sampling the position distribution
Positions = this.initialPositionSampling();
this.InitialPositions = this.initialPositionSampling();
% - sampling the velocity distribution
Velocities = this.initialVelocitySampling();
this.InitialVelocities = this.initialVelocitySampling();
%% Solve ODE
progressbar = Helper.parforNotifications();
progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']);
% calculate the final position of the atoms
DynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
ParticleDynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
Positions = this.InitialPositions;
Velocities = this.InitialVelocities;
parfor Index = 1:this.NumberOfAtoms
DynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
progressbar.PB_iterate();
end
clear Index
this.ParticleDynamicalQuantities = DynamicalQuantities;
%% Calculate the Loading Rate
[LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate();
[LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ParticleDynamicalQuantities);
end

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

@ -11,7 +11,7 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelo
ParticleDynamicalQuantities(i,1:3) = InitialPosition;
ParticleDynamicalQuantities(i,4:6) = InitialVelocity;
rt = InitialPosition;
vt = InitialVelocity;
@ -23,7 +23,7 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelo
ga2 = this.calculateTotalAcceleration(rt,vt) + g;
gv2 = vt .* this.TimeStep;
rt = rt + 0.5 * gv2;
vt = vt + 0.5 * ga2 .* this.TimeStep;
vt = vt + 0.5 *ga2 .* this.TimeStep;
ga3 = this.calculateTotalAcceleration(rt,vt) + g;
gv3 = vt .* this.TimeStep;
@ -35,6 +35,7 @@ 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;
end
end

28
MOT Capture Process Simulation/test_MOTSimulator.m
View File

@ -7,9 +7,9 @@ OptionsStruct.SimulationMode = '2D';
OptionsStruct.TimeStep = 50e-06; % in s
OptionsStruct.SimulationTime = 4e-03; % in s
OptionsStruct.SpontaneousEmission = true;
OptionsStruct.Sideband = false;
OptionsStruct.PushBeam = false;
OptionsStruct.Gravity = false;
OptionsStruct.Sideband = true;
OptionsStruct.PushBeam = true;
OptionsStruct.Gravity = true;
OptionsStruct.BackgroundCollision = true;
OptionsStruct.SaveData = false;
OptionsStruct.SaveDirectory = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation';
@ -43,11 +43,11 @@ clear OptionsStruct
%% - Plot initial distribution
Simulator.setInitialConditions();
% - sampling the position distribution
InitialPositions = Simulator.initialPositionSampling();
Simulator.InitialPositions = Simulator.initialPositionSampling();
% - sampling the velocity distribution
InitialVelocities = Simulator.initialVelocitySampling();
Simulator.InitialVelocities = Simulator.initialVelocitySampling();
NumberOfBins = 100;
Plotting.plotPositionAndVelocitySampling(NumberOfBins, InitialPositions, InitialVelocities);
Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins);
%% - Plot distributions of magnitude and direction of initial velocities
NumberOfBins = 50;
@ -75,28 +75,20 @@ Plotting.plotAngularDistributionForDifferentBeta(Simulator, Beta)
Simulator.setInitialConditions();
Plotting.plotCaptureVelocityVsAngle(Simulator);
%% - Plot Phase Space
%% - Plot Phase Space with Acceleration Field
Simulator.NumberOfAtoms = 100;
MinimumVelocity = 0;
Simulator.NumberOfAtoms = 200;
MaximumVelocity = 150;
NumberOfBins = 200; %Along each axis
IncidentAtomDirection = 0*2*pi/360;
IncidentAtomPosition = 0;
Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
%% - Plot trajectories along the 3 directions
Simulator.NumberOfAtoms = 100;
MaximumVelocity = 150;
IncidentAtomDirection = 0*2*pi/360;
IncidentAtomPosition = 0;
Plotting.plotDynamicalQuantities(Simulator, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition);
%% - Scan parameters
% ONE-PARAMETER SCAN
NumberOfPointsForParam = 5; %iterations of the simulation
NumberOfPointsForParam = 10; %iterations of the simulation
% Scan Cooling Beam Power
PowerArray = linspace(0.1, 1.0, NumberOfPointsForParam) * Simulator.TotalPower;
% Scan Cooling Beam Detuning