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

Newly renamed, restructured Simulator code.

Browse Source
This commit is contained in:
Karthik 2024-06-12 19:09:15 +02:00
parent 41f79671b3
commit 33390f8230
77 changed files with 5705 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
Dipolar Gas Simulator/+Helper/PhysicsConstants.m Normal file
View File

@ -0,0 +1,46 @@
classdef PhysicsConstants < handle
properties (Constant)
% CODATA
PlanckConstant=6.62607015E-34;
PlanckConstantReduced=6.62607015E-34/(2*pi);
FineStructureConstant=7.2973525698E-3;
ElectronMass=9.10938291E-31;
GravitationalConstant=6.67384E-11;
ProtonMass=1.672621777E-27;
AtomicMassUnit=1.66053878283E-27;
BohrRadius=0.52917721092E-10;
BohrMagneton=927.400968E-26;
BoltzmannConstant=1.380649E-23;
StandardGravityAcceleration=9.80665;
SpeedOfLight=299792458;
StefanBoltzmannConstant=5.670373E-8;
ElectronCharge=1.602176634E-19;
VacuumPermeability=1.25663706212E-6;
DielectricConstant=8.8541878128E-12;
ElectronGyromagneticFactor=-2.00231930436153;
AvogadroConstant=6.02214076E23;
ZeroKelvin = 273.15;
GravitationalAcceleration = 9.80553;
% Dy specific constants
Dy164Mass = 163.929174751*1.66053878283E-27;
Dy164IsotopicAbundance = 0.2826;
BlueWavelength = 421.291e-9;
BlueLandegFactor = 1.22;
BlueLifetime = 4.94e-9;
BlueLinewidth = 1/4.94e-9;
RedWavelength = 626.086e-9;
RedLandegFactor = 1.29;
RedLifetime = 1.2e-6;
RedLinewidth = 1/1.2e-6;
PushBeamWaveLength = 626.086e-9;
PushBeamLifetime = 1.2e-6;
PushBeamLinewidth = 1/1.2e-6;
end
methods
function pc = PhysicsConstants()
end
end
end

148
Dipolar Gas Simulator/+Helper/parforNotifications.m Normal file
View File

@ -0,0 +1,148 @@
% Copyright (c) 2019 Andrea Alberti
%
% All rights reserved.
classdef parforNotifications < handle
properties
N; % number of iterations
text = 'Please wait ...'; % text to show
width = 50;
showWarning = true;
end
properties (GetAccess = public, SetAccess = private)
n;
end
properties (Access = private)
inProgress = false;
percent;
DataQueue;
usePercent;
Nstr;
NstrL;
lastComment;
end
methods
function this = parforNotifications()
this.DataQueue = parallel.pool.DataQueue;
afterEach(this.DataQueue, @this.updateStatus);
end
% Start progress bar
function PB_start(this,N,varargin)
assert(isscalar(N) && isnumeric(N) && N == floor(N) && N>0, 'Error: ''N'' must be a scalar positive integer.');
this.N = N;
p = inputParser;
addParameter(p,'message','Please wait: ');
addParameter(p,'usePercentage',true);
parse(p,varargin{:});
this.text = p.Results.message;
assert(ischar(this.text), 'Error: ''Message'' must be a string.');
this.usePercent = p.Results.usePercentage;
assert(isscalar(this.usePercent) && islogical(this.usePercent), 'Error: ''usePercentage'' must be a logical scalar.');
this.percent = 0;
this.n = 0;
this.lastComment = '';
if this.usePercent
fprintf('%s [%s]: %3d%%\n',this.text, char(32*ones(1,this.width)),0);
else
this.Nstr = sprintf('%d',this.N);
this.NstrL = numel(this.Nstr);
fprintf('%s [%s]: %s/%s\n',this.text, char(32*ones(1,this.width)),[char(32*ones(1,this.NstrL-1)),'0'],this.Nstr);
end
this.inProgress = true;
end
% Iterate progress bar
function PB_iterate(this,str)
if nargin == 1
send(this.DataQueue,'');
else
send(this.DataQueue,str);
end
end
function warning(this,warn_id,msg)
if this.showWarning
msg = struct('Action','Warning','Id',warn_id,'Message',msg);
send(this.DataQueue,msg);
end
end
function PB_reprint(this)
p = round(100*this.n/this.N);
this.percent = p;
cursor_pos=1+round((this.width-1)*p/100);
if p < 100
sep_char = '|';
else
sep_char = '.';
end
if this.usePercent
fprintf('%s [%s%s%s]: %3d%%\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p);
else
nstr=sprintf('%d',this.n);
fprintf('%s [%s%s%s]: %s/%s\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr);
end
end
function updateStatus(this,data)
if ischar(data)
this.n = this.n + 1;
p = round(100*this.n/this.N);
if p >= this.percent+1 || this.n == this.N
this.percent = p;
cursor_pos=1+round((this.width-1)*p/100);
if p < 100
sep_char = '|';
else
sep_char = '.';
end
if ~isempty(data)
comment = [' (',data,')'];
else
comment = '';
end
if this.usePercent
fprintf('%s%s%s%s]: %3d%%%s\n',char(8*ones(1,58+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p,comment);
else
nstr=sprintf('%d',this.n);
fprintf('%s%s%s%s]: %s/%s%s\n',char(8*ones(1,55+2*numel(this.Nstr)+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr,comment)
end
this.lastComment = comment;
if p == 100
this.inProgress = false;
end
end
else
switch data.Action
case 'Warning'
warning(data.Id,[data.Message,newline]);
if this.inProgress
this.PB_reprint();
end
end
end
end
end
end

77
Dipolar Gas Simulator/+Plotter/makeMovie.m Normal file
View File

@ -0,0 +1,77 @@
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
RunIdx = 1;
FileDir = dir(sprintf('./Data/Run_%03i/TimeEvolution/*.mat',RunIdx));
NumFiles = numel(FileDir);
QuenchSettings = load(sprintf('./Data/Run_%03i/QuenchSettings',RunIdx),'Quench','Params','Transf','VDk','V');
Transf = QuenchSettings.Transf; Params = QuenchSettings.Params;
x = Transf.x; y = Transf.y; z = Transf.z;
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
mkdir(sprintf('./Data/Run_%03i/Figures',RunIdx))
outputVideo = VideoWriter(fullfile('./Data/Movie.avi'));
outputVideo.FrameRate = 10;
open(outputVideo)
figure(1);
x0 = 800;
y0 = 200;
width = 800;
height = 600;
set(gcf,'position',[x0,y0,width,height])
EVecTemp = [];
for ii = 2:(NumFiles-1)
load(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',RunIdx,ii),'psi','muchem','T','Observ','t_idx');
%Plotting
subplot(2,3,1)
n = abs(psi).^2;
nxz = squeeze(trapz(n*dy,2));
nyz = squeeze(trapz(n*dx,1));
nxy = squeeze(trapz(n*dz,3));
plotxz = pcolor(x,z,nxz'); shading interp
set(plotxz, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
subplot(2,3,2)
plotyz = pcolor(y,z,nyz'); shading interp
set(plotyz, 'EdgeColor', 'none');
xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
subplot(2,3,3)
plotxy = pcolor(x,y,nxy'); shading interp
set(plotxy, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
subplot(2,3,4)
plot(Observ.tVecPlot*1000/Params.w0,Observ.NormVec,'-b')
ylabel('Normalization'); xlabel('$t$ [$m$s]');
subplot(2,3,5)
plot(Observ.tVecPlot*1000/Params.w0,1-2*Observ.PCVec/pi,'-b')
ylabel('Coherence'); xlabel('$t$ [$m$s]');
ylim([0,1])
subplot(2,3,6)
plot(Observ.tVecPlot*1000/Params.w0,Observ.EVec,'-b')
ylabel('E'); xlabel('$t$ [$m$s]');
tVal = Observ.tVecPlot(end)*1000/Params.w0;
sgtitle(sprintf('$\\mu =%.3f \\hbar\\omega_0$, $T=%.1f$nK, $t=%.1f$ms',muchem,T,tVal))
drawnow
saveas(gcf,sprintf('./Data/Run_%03i/Figures/Image_%i.jpg',RunIdx,ii))
img = imread(sprintf('./Data/Run_%03i/Figures/Image_%i.jpg',RunIdx,ii));
writeVideo(outputVideo,img)
% hold off;
clf
end
close(outputVideo)
close(figure(1))
delete(sprintf('./Data/Run_%03i/Figures/*.jpg',RunIdx)) % deleting images after movie is made

47
Dipolar Gas Simulator/+Plotter/plotLive.m Normal file
View File

@ -0,0 +1,47 @@
function LPlot = LivePlot(psi,Params,Transf,Observ)
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
format long
x = Transf.x*Params.l0*1e6;
y = Transf.y*Params.l0*1e6;
z = Transf.z*Params.l0*1e6;
%percentcomplete = linspace(0,1,Params.cut_off/200);
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
%Plotting
subplot(2,3,1)
n = abs(psi).^2;
nxz = squeeze(trapz(n*dy,2));
nyz = squeeze(trapz(n*dx,1));
nxy = squeeze(trapz(n*dz,3));
plotxz = pcolor(x,z,nxz');
set(plotxz, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
subplot(2,3,2)
plotyz = pcolor(y,z,nyz');
set(plotyz, 'EdgeColor', 'none');
xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
subplot(2,3,3)
plotxy = pcolor(x,y,nxy');
set(plotxy, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
subplot(2,3,4)
plot(-log10(Observ.residual),'-b')
ylabel('$-\mathrm{log}_{10}(r)$'); xlabel('steps');
subplot(2,3,5)
plot(Observ.EVec,'-b')
ylabel('$E$'); xlabel('steps');
subplot(2,3,6)
plot(Observ.mucVec,'-b')
ylabel('$\mu$'); xlabel('steps');
% xlim([0,1]); ylim([0,8]);
% xlim([0,1]); ylim([0,8]);

94
Dipolar Gas Simulator/+Plotter/visualizeSpace.m Normal file
View File

@ -0,0 +1,94 @@
function visualizeSpace(Transf)
x = Transf.x;
y = Transf.y;
z = Transf.z;
[X,Y] = meshgrid(x,y);
height = 20;
width = 45;
figure(1)
clf
set(gcf, 'Units', 'centimeters')
set(gcf, 'Position', [2 4 width height])
set(gcf, 'PaperPositionMode', 'auto')
subplot(1,2,1)
zlin = ones(size(X, 1)) * z(1); % Generate z data
mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface
hold on
for idx = 2:length(z)
zlin = ones(size(X, 1))* z(idx); % Generate z data
mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface
end
% set the axes labels' properties
xlabel(gca, {'$x / l_o ~ (m)$'}, ...
'Interpreter', 'latex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
ylabel(gca, {'$y / l_o ~ (m)$'}, ...
'Interpreter', 'latex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
zlabel(gca, {'$z / l_o ~ (m)$'}, ...
'Interpreter', 'latex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
title(gca, 'Real Space', ...
'Interpreter', 'tex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
x = Transf.kx;
y = Transf.ky;
z = Transf.kz;
[X,Y] = meshgrid(x,y);
subplot(1,2,2)
zlin = ones(size(X, 1)) * z(1); % Generate z data
mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface
hold on
for idx = 2:length(z)
zlin = ones(size(X, 1))* z(idx); % Generate z data
mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface
end
% set the axes labels' properties
xlabel(gca, {'$k_x / l_o ~ (m^{-1})$'}, ...
'Interpreter', 'latex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
ylabel(gca, {'$k_y / l_o ~ (m^{-1})$'}, ...
'Interpreter', 'latex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
zlabel(gca, {'$k_z / l_o ~ (m^{-1})$'}, ...
'Interpreter', 'latex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
title(gca, 'Fourier Space', ...
'Interpreter', 'tex', ...
'FontName', 'Times New Roman', ...
'FontSize', 14, ...
'FontWeight', 'normal', ...
'FontAngle', 'normal')
end

50
Dipolar Gas Simulator/+Scripts/analyze.m Normal file
View File

@ -0,0 +1,50 @@
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
format long
runIdx = 6;
load(sprintf('./Data/Run_%03i/psi_gs.mat',runIdx),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
x = Transf.x*Params.l0*1e6;
y = Transf.y*Params.l0*1e6;
z = Transf.z*Params.l0*1e6;
%percentcomplete = linspace(0,1,Params.cut_off/200);
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
%Plotting
subplot(2,3,1)
n = abs(psi).^2;
nxz = squeeze(trapz(n*dy,2));
nyz = squeeze(trapz(n*dx,1));
nxy = squeeze(trapz(n*dz,3));
plotxz = pcolor(x,z,nxz');
set(plotxz, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
subplot(2,3,2)
plotyz = pcolor(y,z,nyz');
set(plotyz, 'EdgeColor', 'none');
xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
subplot(2,3,3)
plotxy = pcolor(x,y,nxy');
set(plotxy, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
subplot(2,3,4)
plot(-log10(Observ.residual),'-b')
ylabel('$-\mathrm{log}_{10}(r)$'); xlabel('steps');
subplot(2,3,5)
plot(Observ.EVec,'-b')
ylabel('$E$'); xlabel('steps');
subplot(2,3,6)
plot(Observ.mucVec,'-b')
ylabel('$\mu$'); xlabel('steps');
% xlim([0,1]); ylim([0,8]);
% xlim([0,1]); ylim([0,8]);
Ecomp = energy_components(psi,Params,Transf,VDk,V);

28
Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m Normal file
View File

@ -0,0 +1,28 @@
function muchem = ChemicalPotential(psi,Params,Transf,VDk,V)
%Parameters
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
% DDIs
frho=fftn(abs(psi).^2);
Phi=real(ifftn(frho.*VDk));
Eddi = (Params.gdd*Phi.*abs(psi).^2);
%Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
%Potential energy
Epot = V.*abs(psi).^2;
%Contact interactions
Eint = Params.gs*abs(psi).^4;
%Quantum fluctuations
Eqf = Params.gammaQF*abs(psi).^5;
%Total energy
muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
muchem = muchem / Params.N;

35
Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m Normal file
View File

@ -0,0 +1,35 @@
function E = EnergyComponents(psi,Params,Transf,VDk,V)
%Parameters
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
% DDIs
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
%Kinetic energy
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
% Potential energy
Epot = V.*abs(psi).^2;
E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
%Contact interactions
Eint = 0.5*Params.gs*abs(psi).^4;
E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
%Quantum fluctuations
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
% plot(Transf.x,abs(psi(:,end/2,end/2+1)).^2)

24
Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m Normal file
View File

@ -0,0 +1,24 @@
function res = NormalizedResiduals(psi,Params,Transf,VDk,V,muchem)
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
% DDIs
frho=fftn(abs(psi).^2);
Phi=real(ifftn(frho.*VDk));
Eddi = Params.gdd*Phi.*psi;
%Kinetic energy
Ekin = ifftn(KEop.*fftn(psi));
%Potential energy
Epot = V.*psi;
%Contact interactions
Eint = Params.gs*abs(psi).^2.*psi;
%Quantum fluctuations
Eqf = Params.gammaQF*abs(psi).^3.*psi;
%Total energy
res = trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);

39
Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m Normal file
View File

@ -0,0 +1,39 @@
function VDkSemi = NumericalHankelTransform(kr, kz, Rmax, Zmax, Nr)
% accuracy inputs for numerical integration
if(nargin==4)
Nr = 5e4;
end
Nz = 64;
farRmultiple = 2000;
% midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
dr=(farRmultiple-1)*Rmax/Nr;
r = ((1:Nr)'-0.5)*dr+Rmax;
dz=Zmax/Nz;
z = ((1:Nz)-0.5)*dz;
[R, Z] = ndgrid(r,z);
Rsq = R.^2 + Z.^2;
% real space interaction to be transformed
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
% do the Hankel/Fourier-Bessel transform numerically
% prestore to ensure each besselj is only calculated once
% cell is faster than (:,:,krn) slicing
Nkr = numel(kr);
besselr = cell(Nkr,1);
for krn = 1:Nkr
besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
end
VDkSemi = zeros([Nkr,numel(kz)]);
for kzn = 1:numel(kz)
igrandbasez = repmat(cos(kz(kzn)*z),Nr,1) .* igrandbase;
for krn = 1:Nkr
igrand = igrandbasez.*besselr{krn};
VDkSemi(krn,kzn) = VDkSemi(krn,kzn) - sum(igrand(:))*dz*dr;
end
end

58
Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m Normal file
View File

@ -0,0 +1,58 @@
function [m_Order] = OrderParameter(psi,Transf,Params,VDk,V,T,muchem)
NumRealiz = 100;
Mx = numel(Transf.x);
My = numel(Transf.y);
Mz = numel(Transf.z);
r = normrnd(0,1,size(psi));
theta = rand(size(psi));
noise = r.*exp(2*pi*1i*theta);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Gamma = 1-1i*Params.gamma_S;
dt = Params.dt;
avgpsi = 0;
avgpsi2 = 0;
for jj = 1:NumRealiz
%generate initial state
xi = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
swapx = randi(length(Transf.x),1,length(Transf.x));
swapy = randi(length(Transf.y),1,length(Transf.y));
swapz = randi(length(Transf.z),1,length(Transf.z));
psi_j = psi + xi * noise(swapx,swapy,swapz);
% --- % propagate forward in time 1 time step:
%kin
psi_j = fftn(psi_j);
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
psi_j = ifftn(psi_j);
%DDI
frho = fftn(abs(psi_j).^2);
Phi = real(ifftn(frho.*VDk));
%Real-space
psi_j = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));
%kin
psi_j = fftn(psi_j);
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
psi_j = ifftn(psi_j);
%Projection
kcut = sqrt(2*Params.e_cut);
K = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)<kcut.^2;
psi_j = ifftn(K.*fftn(psi_j));
% --- %
avgpsi = avgpsi + abs(sum(psi_j(:)))/NumRealiz;
avgpsi2 = avgpsi2 + sum(abs(psi_j(:)).^2)/NumRealiz;
end
m_Order = 1/sqrt(Mx*My*Mz)*avgpsi/sqrt(avgpsi2);
end

18
Dipolar Gas Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m Normal file
View File

@ -0,0 +1,18 @@
function [PhaseC] = PhaseCoherence(psi,Transf,Params)
norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
psi = psi/sqrt(norm);
NumGlobalShifts = 800;
betas = []; phishift = [];
for jj = 1:NumGlobalShifts
phishift(jj) = -pi + 2*pi*(jj-1)/NumGlobalShifts;
betas(jj) = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
end
[minbeta,minidx] = min(betas);
psi = psi*exp(-1i*phishift(minidx));
phi = angle(psi);
avgphi = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
PhaseC = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;

31
Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m Normal file
View File

@ -0,0 +1,31 @@
function E = TotalEnergy(psi,Params,Transf,VDk,V)
%Parameters
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
% DDIs
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
%Kinetic energy
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
% Potential energy
Epot = V.*abs(psi).^2;
%Contact interactions
Eint = 0.5*Params.gs*abs(psi).^4;
%Quantum fluctuations
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
E = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;

33
Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m Normal file
View File

@ -0,0 +1,33 @@
function VDk = VDCutoff(Params, Transf)
% makes the dipolar interaction matrix, size numel(Params.kr) * numel(Params.kz)
% Rmax and Zmax are the interaction cutoffs
% VDk needs to be multiplied by Cdd
% approach is that of Lu, PRA 82, 023622 (2010)
Zcutoff = Params.Lz/2;
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
% Analytic part of cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
cossq = cos(alph).^2;
VDk = cossq-1/3;
sinsq = 1 - cossq;
VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
% Nonanalytic part
% For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
Params.Lr = 0.5*min(Params.Lx,Params.Ly);
Params.Nr = max(Params.Nx,Params.Ny);
[TransfRad] = SetupSpaceRadial(Params);
VDkNon = NumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
% Interpolating the nonanalytic part onto 3D grid
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
KR3D = sqrt(KX3D.^2 + KY3D.^2);
fullVDK = [flip(VDkNon',2),VDkNon']';
VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
VDkNon = fftshift(VDkNon);
VDk = VDk + VDkNon;

22
Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m Normal file
View File

@ -0,0 +1,22 @@
function [psi,V,VDk] = Initialize(Params,Transf)
format long
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
% == Trap potential == %
V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2);
% == Calculating the DDIs == %
if isfile('./Data/VDk_M.mat')
VDk = load(sprintf('./Data/VDk_M.mat'));
VDk = VDk.VDk;
else
VDk = VDCutoff(Params, Transf);
save(sprintf('./Data/VDk_M.mat'),'VDk');
end
disp('Finished DDI')
% == Setting up the initial wavefunction == %
psi = SetupWavefunction();
end

35
Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m Normal file
View File

@ -0,0 +1,35 @@
%-% Run Simulation %-%
clearvars
cd(fileparts(matlab.desktop.editor.getActiveFilename));
pwd
% --- Obtain simulation parameters ---
[Params] = SetupParameters();
% --- Set up spatial grids and transforms ---
[Transf] = SetupSpace(Params);
% Plotter.visualizeSpace(Transf)
%%
% --- Initialize ---
mkdir(sprintf('./Data'))
[psi,V,VDk] = Initialize(Params,Transf);
Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
t_idx = 1; %Start at t = 0;
Observ.res_idx = 1;
%%
% --- Job Settings ---
njob = 6;
mkdir(sprintf('./Data/Run_%03i',njob))
% --- Run Simulation ---
% Imaginary Time Evolution
[psi] = SplitStepFourierImaginaryTime(psi,Params,Transf,VDk,V,njob,t_idx,Observ);

101
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m Normal file
View File

@ -0,0 +1,101 @@
function [Params] = SetupParameters()
%%--%% Parameters %%--%%
%========= Simulation =========%
pert = 0; % 0 = no perturbation during real-time, 1=perturbation
%method=1; % 0 = normal dipolar potential, 1=spherical cut-off, 2=cylindrical cut-off
% Tolerances
Params.Etol = 5e-10;
Params.rtol = 1e-5;
Params.cut_off = 2e6; % sometimes the imaginary time gets a little stuck
% even though the solution is good, this just stops it going on forever
%========= Constants =========%
hbar = 1.0545718e-34; % Planck constant [J.s]
kbol = 1.38064852e-23; % Boltzmann Constant [J/K]
mu0 = 1.25663706212e-6; % Vacuum Permeability [N/A^2] --
muB = 9.274009994e-24; % Bohr Magneton [J/T]
a0 = 5.2917721067e-11; % Bohr radius [m]
m0 = 1.660539066e-27; % Atomic mass [kg]
w0 = 2*pi*100; % Angular frequency unit [s^-1]
mu0factor = 0.3049584233607396; % =(m0/me)*pi*alpha^2 -- me=mass of electron, alpha=fine struct. const.
% mu0=mu0factor *hbar^2*a0/(m0*muB^2)
%=============================%
% Number of points in each direction
Params.Nx = 64;
Params.Ny = 64;
Params.Nz = 48;
% Dimensions (in units of l0)
Params.Lx = 40;
Params.Ly = 40;
Params.Lz = 20;
% Masses
Params.m = 164*m0;
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
% Atom numbers
% Params.ppum = 2500; % particles per micron
% Params.N = Params.Lz*Params.ppum*l0*1e6;
Params.N = 10^6;
% Dipole angle
Params.theta = pi/2; % pi/2 dipoles along x, theta=0 dipoles along z
% Dipole lengths (units of muB)
Params.mu = 9.93*muB;
% Scattering lengths
Params.as = 86*a0;
% Trapping frequencies
Params.wx = 2*pi*125;
Params.wy = 2*pi*125;
Params.wz = 2*pi*250;
% Time step
Params.dt = 0.0005;
Params.mindt = 1e-6; %Minimum size for a time step using adaptive dt
% Stochastic GPE
Params.gamma_S = 7.5*10^(-3); % gamma for the stochastic GPE
Params.muchem = 12.64*Params.wz/w0; % fixing the chemical potential for the stochastic GPE
% ================ Parameters defined by those above ================ %
% Contact interaction strength (units of l0/m)
Params.gs = 4*pi*Params.as/l0;
% Dipole lengths
Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2);
% DDI strength
Params.gdd = 12*pi*Params.add/l0; %sometimes the 12 is a 4? --> depends on how Vdk (DDI) is defined
% Trap gamma
Params.gx=(Params.wx/w0)^2;
Params.gy=(Params.wy/w0)^2;
Params.gz=(Params.wz/w0)^2;
% == Calculating quantum fluctuations == %
eps_dd = Params.add/Params.as;
if eps_dd == 0
Q5 = 1;
elseif eps_dd == 1
Q5 = 3*sqrt(3)/2;
else
yeps = (1-eps_dd)/(3*eps_dd);
Q5 = (3*eps_dd)^(5/2)*( (8+26*yeps+33*yeps^2)*sqrt(1+yeps) + 15*yeps^3*log((1+sqrt(1+yeps))/sqrt(yeps)) )/48;
Q5 = real(Q5);
end
Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5;
% Loading the rest into Params
Params.hbar = hbar; Params.kbol = kbol; Params.mu0 = mu0; Params.muB = muB; Params.a0 = a0;
Params.w0 = w0; Params.l0 = l0;
end

33
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m Normal file
View File

@ -0,0 +1,33 @@
function [Transf] = SetupSpace(Params)
Transf.Xmax = 0.5*Params.Lx;
Transf.Ymax = 0.5*Params.Ly;
Transf.Zmax = 0.5*Params.Lz;
Nz = Params.Nz; Nx = Params.Nx; Ny = Params.Ny;
% Fourier grids
x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
Kmax = pi*Params.Nx/Params.Lx;
kx = linspace(-Kmax,Kmax,Nx+1);
kx = kx(1:end-1); dkx = kx(2)-kx(1);
kx = fftshift(kx);
y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
Kmax = pi*Params.Ny/Params.Ly;
ky = linspace(-Kmax,Kmax,Ny+1);
ky = ky(1:end-1); dky = ky(2)-ky(1);
ky = fftshift(ky);
z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz);
Kmax = pi*Params.Nz/Params.Lz;
kz = linspace(-Kmax,Kmax,Nz+1);
kz = kz(1:end-1); dkz = kz(2)-kz(1);
kz = fftshift(kz);
[Transf.X,Transf.Y,Transf.Z]=ndgrid(x,y,z);
[Transf.KX,Transf.KY,Transf.KZ]=ndgrid(kx,ky,kz);
Transf.x = x; Transf.y = y; Transf.z = z;
Transf.kx = kx; Transf.ky = ky; Transf.kz = kz;
Transf.dx = x(2)-x(1); Transf.dy = y(2)-y(1); Transf.dz = z(2)-z(1);
Transf.dkx = dkx; Transf.dky = dky; Transf.dkz = dkz;
end

307
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m Normal file
View File

@ -0,0 +1,307 @@
function [Transf] = SetupSpaceRadial(Params,morder)
Zmax = 0.5*Params.Lz;
Rmax = Params.Lr;
Nz = Params.Nz;
Nr = Params.Nr;
if(nargin==1)
morder=0; %only do Bessel J0
end
% Fourier grids
z=linspace(-Zmax,Zmax,Nz+1);
z=z(1:end-1);
dz=z(2)-z(1);
Kmax=Nz*2*pi/(4*Zmax);
kz=linspace(-Kmax,Kmax,Nz+1);
kz=kz(1:end-1);
% Hankel grids and transform
H = hankelmatrix(morder,Rmax,Nr);
r=H.r(:);
kr=H.kr(:);
T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi);
Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi);
wr=H.wr;
wk=H.wk;
% H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J)
[Transf.R,Transf.Z]=ndgrid(r,z);
[Transf.KR,Transf.KZ]=ndgrid(kr,kz);
Transf.T=T;
Transf.Tinv=Tinv;
Transf.r=r;
Transf.kr=kr;
Transf.z=z;
Transf.kz=kz;
Transf.wr=wr;
Transf.wk=wk;
Transf.Rmax=Rmax;
Transf.Zmax=Zmax;
Transf.dz=z(2)-z(1);
Transf.dkz=kz(2)-kz(1);
%b1=Transf;
function s_HT = hankelmatrix(order, rmax, Nr, eps_roots)
%HANKEL_MATRIX: Generates data to use for Hankel Transforms
%
% s_HT = hankel_matrix(order, rmax, Nr, eps_roots)
%
% s_HT = Structure containing data to use for the pQDHT
% order = Transform order
% rmax = Radial extent of transform
% Nr = Number of sample points
% eps_roots = Error in estimation of roots of Bessel function (optional)
%
% s_HT:
% order, rmax, Nr = As above
% J_roots = Roots of the pth order Bessel fn.
% J_roots_N1 = (N+1)th root
% r = Radial co-ordinate vector
% v = frequency co-ordinate vector
% kr = Radial wave number co-ordinate vector
% vmax = Limiting frequency
% = roots_N1 / (2*pi*rmax)
% S = rmax * 2*pi*vmax (S product)
% T = Transform matrix
% J = Scaling vector
% = J_(order+1){roots}
%
% The algorithm used is that from:
% "Computation of quasi-discrete Hankel transforms of the integer
% order for propagating optical wave fields"
% Manuel Guizar-Sicairos and Julio C. Guitierrez-Vega
% J. Opt. Soc. Am. A 21(1) 53-58 (2004)
%
% The algorithm also calls the function:
% zn = bessel_zeros(1, p, Nr+1, 1e-6),
% where p and N are defined above, to calculate the roots of the bessel
% function. This algorithm is taken from:
% "An Algorithm with ALGOL 60 Program for the Computation of the
% zeros of the Ordinary Bessel Functions and those of their
% Derivatives".
% N. M. Temme
% Journal of Computational Physics, 32, 270-279 (1979)
%
% Example: Propagation of radial field
%
% % Note the use of matrix and element products / divisions
% H = hankel_matrix(0, 1e-3, 512);
% DR0 = 50e-6;
% Ur0 = exp(-(H.r/DR0).^2);
% Ukr0 = H.T * (Ur0./H.J);
% k0 = 2*pi/800e-9;
% kz = realsqrt((k0^2 - H.kr.^2).*(k0>H.kr));
% z = (-5e-3:1e-5:5e-3);
% Ukrz = (Ukr0*ones(1,length(z))).*exp(i*kz*z);
% Urz = (H.T * Ukrz) .* (H.J * ones(1,length(z)));
%
% See also bessel_zeros, besselj
if (~exist('eps_roots', 'var')||isemtpy(eps_roots))
s_HT.eps_roots = 1e-6;
else
s_HT.eps_roots = eps_roots;
end
s_HT.order = order;
s_HT.rmax = rmax;
s_HT.Nr = Nr;
% Calculate N+1 roots:
J_roots = bessel_zeros(1, s_HT.order, s_HT.Nr+1, s_HT.eps_roots);
s_HT.J_roots = J_roots(1:end-1);
s_HT.J_roots_N1 = J_roots(end);
% Calculate co-ordinate vectors
s_HT.r = s_HT.J_roots * s_HT.rmax / s_HT.J_roots_N1;
s_HT.v = s_HT.J_roots / (2*pi * s_HT.rmax);
s_HT.kr = 2*pi * s_HT.v;
s_HT.kmax = s_HT.J_roots_N1 / (s_HT.rmax);
s_HT.vmax = s_HT.J_roots_N1 / (2*pi * s_HT.rmax);
s_HT.S = s_HT.J_roots_N1;
% Calculate hankel matrix and vectors
% I use (p=order) and (p1=order+1)
Jp = besselj(s_HT.order, (s_HT.J_roots) * (s_HT.J_roots.') / s_HT.S);
Jp1 = abs(besselj(s_HT.order+1, s_HT.J_roots));
s_HT.T = 2*Jp./(Jp1 * (Jp1.') * s_HT.S);
s_HT.J = Jp1;
s_HT.wr=2./((s_HT.kmax)^2*abs(Jp1).^2);
s_HT.wk=2./((s_HT.rmax)^2*abs(Jp1).^2);
return
function z = bessel_zeros(d, a, n, e)
%BESSEL_ZEROS: Finds the first n zeros of a bessel function
%
% z = bessel_zeros(d, a, n, e)
%
% z = zeros of the bessel function
% d = Bessel function type:
% 1: Ja
% 2: Ya
% 3: Ja'
% 4: Ya'
% a = Bessel order (a>=0)
% n = Number of zeros to find
% e = Relative error in root
%
% This function uses the routine described in:
% "An Algorithm with ALGOL 60 Program for the Computation of the
% zeros of the Ordinary Bessel Functions and those of their
% Derivatives".
% N. M. Temme
% Journal of Computational Physics, 32, 270-279 (1979)
z = zeros(n, 1);
aa = a^2;
mu = 4*aa;
mu2 = mu^2;
mu3 = mu^3;
mu4 = mu^4;
if (d<3)
p = 7*mu - 31;
p0 = mu - 1;
if ((1+p)==p)
p1 = 0;
q1 = 0;
else
p1 = 4*(253*mu2 - 3722*mu+17869)*p0/(15*p);
q1 = 1.6*(83*mu2 - 982*mu + 3779)/p;
end
else
p = 7*mu2 + 82*mu - 9;
p0 = mu + 3;
if ((p+1)==1)
p1 = 0;
q1 = 0;
else
p1 = (4048*mu4 + 131264*mu3 - 221984*mu2 - 417600*mu + 1012176)/(60*p);
q1 = 1.6*(83*mu3 + 2075*mu2 - 3039*mu + 3537)/p;
end
end
if (d==1)|(d==4)
t = .25;
else
t = .75;
end
tt = 4*t;
if (d<3)
pp1 = 5/48;
qq1 = -5/36;
else
pp1 = -7/48;
qq1 = 35/288;
end
y = .375*pi;
if (a>=3)
bb = a^(-2/3);
else
bb = 1;
end
a1 = 3*a - 8;
% psi = (.5*a + .25)*pi;
for s=1:n
if ((a==0)&(s==1)&(d==3))
x = 0;
j = 0;
else
if (s>=a1)
b = (s + .5*a - t)*pi;
c = .015625/(b^2);
x = b - .125*(p0 - p1*c)/(b*(1 - q1*c));
else
if (s==1)
switch (d)
case (1)
x = -2.33811;
case (2)
x = -1.17371;
case (3)
x = -1.01879;
otherwise
x = -2.29444;
end
else
x = y*(4*s - tt);
v = x^(-2);
x = -x^(2/3) * (1 + v*(pp1 + qq1*v));
end
u = x*bb;
v = fi(2/3 * (-u)^1.5);
w = 1/cos(v);
xx = 1 - w^2;
c = sqrt(u/xx);
if (d<3)
x = w*(a + c*(-5/u - c*(6 - 10/xx))/(48*a*u));
else
x = w*(a + c*(7/u + c*(18 - 14/xx))/(48*a*u));
end
end
j = 0;
while ((j==0)|((j<5)&(abs(w/x)>e)))
xx = x^2;
x4 = x^4;
a2 = aa - xx;
r0 = bessr(d, a, x);
j = j+1;
if (d<3)
u = r0;
w = 6*x*(2*a + 1);
p = (1 - 4*a2)/w;
q = (4*(xx-mu) - 2 - 12*a)/w;
else
u = -xx*r0/a2;
v = 2*x*a2/(3*(aa+xx));
w = 64*a2^3;
q = 2*v*(1 + mu2 + 32*mu*xx + 48*x4)/w;
p = v*(1 + (40*mu*xx + 48*x4 - mu2)/w);
end
w = u*(1 + p*r0)/(1 + q*r0);
x = x+w;
end
z(s) = x;
end
end
function FI = fi(y)
c1 = 1.570796;
if (~y)
FI = 0;
elseif (y>1e5)
FI = c1;
else
if (y<1)
p = (3*y)^(1/3);
pp = p^2;
p = p*(1 + pp*(pp*(27 - 2*pp) - 210)/1575);
else
p = 1/(y + c1);
pp = p^2;
p = c1 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328))))/3465);
end
pp = (y+p)^2;
r = (p - atan(p+y))/pp;
FI = p - (1+pp)*r*(1 + r/(p+y));
end
return
function Jr = bessr(d, a, x)
switch (d)
case (1)
Jr = besselj(a, x)./besselj(a+1, x);
case (2)
Jr = bessely(a, x)./bessely(a+1, x);
case (3)
Jr = a./x - besselj(a+1, x)./besselj(a, x);
otherwise
Jr = a./x - bessely(a+1, x)./bessely(a, x);
end
return

23
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m Normal file
View File

@ -0,0 +1,23 @@
function [psi] = SetupWavefunction()
ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0;
elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0;
Rx = 8*ellx; Ry = 8*elly; Rz = 8*ellz;
X0 = 0.0*Transf.Xmax; Y0 = 0.0*Transf.Ymax; Z0 = 0*Transf.Zmax;
psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2-(Z-Z0).^2/Rz^2);
cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = psi/sqrt(cur_norm);
%add some noise
r = normrnd(0,1,size(X));
theta = rand(size(X));
noise = r.*exp(2*pi*1i*theta);
psi = 0*psi + 1*noise;
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
end

90
Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m Normal file
View File

@ -0,0 +1,90 @@
function [psi] = SplitStepFourier(psi,Params,Transf,VDk,V,njob,t_idx,Observ)
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
dt=-1j*abs(Params.dt);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Observ.residual = 1; Observ.res = 1;
muchem = Simulator.ChemicalPotential(psi,Params,Transf,VDk,V);
AdaptIdx = 0;
while t_idx < Params.cut_off
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%DDI
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
%Real-space
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%Renorm
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
muchem = Simulator.ChemicalPotential(psi,Params,Transf,VDk,V);
if mod(t_idx,1000) == 0
%Change in Energy
E = Simulator.TotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
%Chemical potential
Observ.mucVec = [Observ.mucVec muchem];
%Normalized residuals
res = Simulator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
Observ.residual = [Observ.residual res];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
%Adaptive time step -- Careful, this can quickly get out of control
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
if relres <1e-5
if AdaptIdx > 4 && abs(dt) > Params.mindt
dt = dt / 2;
fprintf('Time step changed to '); disp(dt);
AdaptIdx = 0;
elseif AdaptIdx > 4 && abs(dt) < Params.mindt
break
else
AdaptIdx = AdaptIdx + 1;
end
else
AdaptIdx = 0;
end
end
if any(isnan(psi(:)))
disp('NaNs encountered!')
break
end
t_idx=t_idx+1;
end
%Change in Energy
E = Simulator.TotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Phase coherence
[PhaseC] = Simulator.PhaseCoherence(psi,Transf,Params);
Observ.PCVec = [Observ.PCVec PhaseC];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
end

BIN
MOT Simulator/+Helper/ImageSelection.class Normal file
View File

Binary file not shown.

38
MOT Simulator/+Helper/ImageSelection.java Normal file
View File

@ -0,0 +1,38 @@
/*
* Based on code snippet from
* http://java.sun.com/developer/technicalArticles/releases/data/
*
* Copyright © 2008, 2010 Oracle and/or its affiliates. All rights reserved. Use is subject to license terms.
*/
import java.awt.image.BufferedImage;
import java.awt.datatransfer.*;
public class ImageSelection implements Transferable {
private static final DataFlavor flavors[] =
{DataFlavor.imageFlavor};
private BufferedImage image;
public ImageSelection(BufferedImage image) {
this.image = image;
}
// Transferable
public Object getTransferData(DataFlavor flavor) throws UnsupportedFlavorException {
if (flavor.equals(flavors[0]) == false) {
throw new UnsupportedFlavorException(flavor);
}
return image;
}
public DataFlavor[] getTransferDataFlavors() {
return flavors;
}
public boolean isDataFlavorSupported(DataFlavor
flavor) {
return flavor.equals(flavors[0]);
}
}

46
MOT Simulator/+Helper/PhysicsConstants.m Normal file
View File

@ -0,0 +1,46 @@
classdef PhysicsConstants < handle
properties (Constant)
% CODATA
PlanckConstant=6.62607015E-34;
PlanckConstantReduced=6.62607015E-34/(2*pi);
FineStructureConstant=7.2973525698E-3;
ElectronMass=9.10938291E-31;
GravitationalConstant=6.67384E-11;
ProtonMass=1.672621777E-27;
AtomicMassUnit=1.66053878283E-27;
BohrRadius=0.52917721092E-10;
BohrMagneton=927.400968E-26;
BoltzmannConstant=1.380649E-23;
StandardGravityAcceleration=9.80665;
SpeedOfLight=299792458;
StefanBoltzmannConstant=5.670373E-8;
ElectronCharge=1.602176634E-19;
VacuumPermeability=1.25663706212E-6;
DielectricConstant=8.8541878128E-12;
ElectronGyromagneticFactor=-2.00231930436153;
AvogadroConstant=6.02214076E23;
ZeroKelvin = 273.15;
GravitationalAcceleration = 9.80553;
% Dy specific constants
Dy164Mass = 163.929174751*1.66053878283E-27;
Dy164IsotopicAbundance = 0.2826;
BlueWavelength = 421.291e-9;
BlueLandegFactor = 1.22;
BlueLifetime = 4.94e-9;
BlueLinewidth = 1/4.94e-9;
RedWavelength = 626.086e-9;
RedLandegFactor = 1.29;
RedLifetime = 1.2e-6;
RedLinewidth = 1/1.2e-6;
PushBeamWaveLength = 626.086e-9;
PushBeamLifetime = 1.2e-6;
PushBeamLinewidth = 1/1.2e-6;
end
methods
function pc = PhysicsConstants()
end
end
end

15
MOT Simulator/+Helper/bringFiguresWithTagInForeground.m Normal file
View File

@ -0,0 +1,15 @@
function output = bringFiguresWithTagInForeground()
figure_handles = findobj('type','figure');
for idx = 1:length(figure_handles)
if ~isempty(figure_handles(idx).Tag)
figure(figure_handles(idx));
end
end
if nargout > 0
output = figure_handles;
end
end

10
MOT Simulator/+Helper/calculateDistanceFromPointToLine.m Normal file
View File

@ -0,0 +1,10 @@
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);
%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

6
MOT Simulator/+Helper/convertstruct2cell.m Normal file
View File

@ -0,0 +1,6 @@
function CellOut = convertstruct2cell(StructIn)
% CellOut = Convertstruct2cell(StructIn)
% converts a struct into a cell-matrix where the first column contains
% the fieldnames and the second the contents
CellOut = [fieldnames(StructIn) struct2cell(StructIn)]';
end

18
MOT Simulator/+Helper/findAllZeroCrossings.m Normal file
View File

@ -0,0 +1,18 @@
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

191
MOT Simulator/+Helper/getFigureByTag.m Normal file
View File

@ -0,0 +1,191 @@
function figure_handle = getFigureByTag(tag_name, varargin)
% figure_handle = getFigureByTag(tag_name, varargin)
%
% Example code:
% f_h = getFigureByTag('survivalMeasurement','Name','Survival')
%
% clf(f_h);
% a_h = gca(f_h);
% xlim(a_h,[10,100]);
% % custom position
% f_h.Position = [4052.3 719.67 560 420];
assert(nargin>=1 && ischar(tag_name),'You must specify ``tag_name'' as a string.');
f_h = findobj('type','figure','tag',tag_name);
if isempty(f_h)
f_h = figure('Tag',tag_name,varargin{:});
defaultNewFigProperties = {'Color','w','NumberTitle','off','Name',sprintf('Fig. %d',f_h.Number)};
varargin = [defaultNewFigProperties,varargin];
else
f_h = f_h(1);
end
if ~isempty(varargin)
set(f_h,varargin{:});
end
addCopyButton(f_h);
if nargout > 0
figure_handle = f_h;
else
set(groot,'CurrentFigure',f_h);
end
end
function addCopyButton(f_h)
if(strcmp(f_h.ToolBar,'none'))
return
end
tb = findall(f_h,'Type','uitoolbar');
pt = findall(tb, 'tag', 'Custom.CopyPlot' );
if isempty(pt)
pt = uipushtool(tb);
else
pt = pt(1);
end
cdata = zeros(16,16,3);
% Evernote Logo
% cdata(:,:,1) =[255 NaN NaN NaN NaN 99 11 27 175 NaN NaN NaN NaN NaN NaN 255
% NaN NaN NaN 251 93 14 0 0 0 66 70 106 210 NaN NaN NaN
% NaN NaN NaN 42 0 43 0 0 0 0 0 0 20 185 NaN NaN
% NaN 243 56 0 42 82 0 0 0 0 0 0 0 45 NaN NaN
% NaN 156 44 64 113 65 0 0 0 0 0 0 0 32 NaN NaN
% 136 9 26 28 11 0 0 0 0 0 0 0 0 10 188 NaN
% 132 0 0 0 0 0 0 0 0 0 136 175 16 0 133 NaN
% NaN 28 0 0 0 0 0 0 0 0 152 238 50 0 124 NaN
% NaN 58 0 0 0 0 0 0 0 0 0 9 0 0 71 NaN
% NaN 175 0 0 0 0 0 61 15 0 0 0 0 0 100 NaN
% NaN NaN 143 12 0 0 0 210 195 87 17 0 0 0 126 NaN
% NaN NaN NaN 183 118 50 150 NaN NaN 110 219 78 0 0 160 NaN
% NaN NaN NaN NaN NaN NaN NaN 191 0 35 NaN 150 0 23 NaN NaN
% NaN NaN NaN NaN NaN NaN NaN 124 0 172 NaN 81 0 93 NaN NaN
% 255 NaN NaN NaN NaN NaN NaN 183 0 0 0 0 51 228 NaN 245
% 253 254 NaN NaN NaN NaN NaN NaN 156 63 45 100 NaN NaN 255 255]/255.;
%
%
% cdata(:,:,2) = [255 255 255 255 255 216 166 171 225 229 218 229 247 255 255 255
% 255 255 255 255 201 166 159 157 167 188 189 200 243 255 255 255
% 237 238 255 181 159 183 164 170 163 158 160 157 169 233 248 250
% 224 235 188 140 182 195 161 168 168 168 168 169 147 186 244 240
% 255 226 175 185 207 189 161 168 168 168 168 168 159 179 249 249
% 227 172 172 179 172 163 169 168 168 170 163 155 160 173 231 237
% 215 161 163 165 166 168 168 168 168 162 215 228 172 163 209 219
% 248 178 159 168 168 168 168 168 168 159 220 249 185 158 208 222
% 249 192 151 169 168 168 169 160 163 172 163 159 166 167 194 204
% 246 229 155 157 168 169 159 188 174 154 162 167 166 166 202 214
% 212 231 218 168 157 153 165 255 242 190 171 159 167 166 207 220
% 218 203 251 243 206 181 230 210 208 207 242 196 154 168 223 232
% 255 224 232 250 237 214 244 194 152 178 255 223 145 175 250 252
% 255 255 244 239 222 213 240 214 149 228 254 199 136 203 244 232
% 255 255 255 246 231 246 246 232 165 159 167 147 184 253 254 242
% 253 254 255 255 254 255 255 255 231 183 178 199 249 255 255 255]/255.;
%
%
% cdata(:,:,3) = [255 255 255 255 255 117 38 50 187 211 170 190 234 255 255 255
% 255 254 255 255 120 51 27 20 39 97 98 122 220 255 255 255
% 238 252 246 73 22 71 37 49 35 20 24 18 49 196 231 231
% 232 242 86 0 78 108 29 45 45 45 45 46 0 82 214 201
% 255 175 63 85 139 98 27 45 45 45 45 45 23 72 233 231
% 167 51 57 72 55 32 47 45 45 50 34 14 27 57 201 218
% 154 30 33 38 39 45 45 45 45 31 157 188 53 34 153 180
% 234 67 24 45 45 45 45 44 45 24 169 241 83 20 146 182
% 241 99 4 48 45 45 47 28 35 53 32 26 39 44 104 127
% 238 192 14 20 45 47 27 97 56 10 29 44 41 40 127 158
% 214 253 169 37 20 16 34 218 207 105 55 23 42 40 147 182
% 218 214 241 201 138 71 177 225 181 130 224 107 12 45 175 197
% 255 233 202 218 212 132 230 196 27 61 255 172 0 64 240 242
% 255 255 219 197 176 160 237 143 0 195 245 110 0 123 230 230
% 255 255 255 227 197 241 244 202 36 24 39 0 81 228 242 245
% 253 254 255 255 254 255 255 255 191 78 71 121 221 255 255 255]/255.;
%OneNote logo
cdata(:,:,1) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 245 213 213 213 213 213 213 213 184 184 215 255
255 255 255 255 241 213 213 213 213 213 213 213 184 184 208 255
255 233 204 204 194 176 176 185 213 213 213 213 184 184 208 255
255 154 101 101 101 101 101 103 213 213 213 206 162 162 193 255
255 152 101 183 116 152 115 101 213 213 213 206 162 162 193 255
255 152 101 207 189 178 122 101 213 213 213 206 162 162 193 255
255 152 101 199 152 224 122 101 213 213 213 195 128 128 170 255
255 152 101 166 101 183 115 101 213 213 213 195 128 128 170 255
255 154 101 101 101 101 101 103 213 213 213 195 128 128 170 255
255 233 204 204 194 176 176 185 213 213 213 183 95 95 148 255
255 255 255 255 241 213 213 213 213 213 213 183 94 94 148 255
255 255 255 255 245 213 213 213 213 213 213 183 94 94 163 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.;
cdata(:,:,2) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 219 112 110 110 110 110 110 134 84 84 158 255
255 255 255 255 207 110 110 110 110 110 110 134 84 84 141 255
255 222 178 178 146 81 81 88 110 110 110 134 84 84 141 255
255 102 23 23 23 23 23 24 110 110 110 125 58 58 123 255
255 100 23 147 46 100 44 23 110 110 110 125 58 58 123 255
255 100 23 183 156 139 55 23 110 110 110 125 58 58 123 255
255 100 23 170 99 208 55 23 110 110 110 119 38 38 109 255
255 100 23 121 23 146 44 23 110 110 110 119 38 38 109 255
255 102 23 23 23 23 23 24 110 110 110 119 38 38 109 255
255 222 178 178 146 81 81 88 110 110 110 118 37 37 109 255
255 255 255 255 207 110 110 110 110 110 110 118 37 37 110 255
255 255 255 255 219 112 110 110 110 110 110 118 37 37 131 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.;
cdata(:,:,3) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255 255 255 255 246 229 229 240 255
255 255 255 255 255 255 255 255 255 255 255 246 229 229 238 255
255 242 224 224 224 224 224 232 255 255 255 246 229 229 238 255
255 194 163 163 163 163 163 164 255 255 255 244 223 223 234 255
255 194 163 212 172 194 171 163 255 255 255 244 223 223 234 255
255 194 163 226 216 209 176 163 255 255 255 244 223 223 234 255
255 194 163 221 193 236 176 163 255 255 255 240 209 209 224 255
255 194 163 202 163 212 171 163 255 255 255 240 209 209 224 255
255 194 163 163 163 163 163 164 255 255 255 240 209 209 224 255
255 242 224 224 224 224 224 232 255 255 255 223 161 161 192 255
255 255 255 255 255 255 255 255 255 255 255 223 160 160 192 255
255 255 255 255 255 255 255 255 255 255 255 223 160 160 201 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.;
pt.Tag = 'Custom.CopyPlot';
pt.CData = cdata;
pt.Separator = true;
pt.ClickedCallback = @copyToClipboard;
end
function copyToClipboard(~,~)
fig_h = get(get(gcbo,'Parent'),'Parent');
if strcmp(fig_h.WindowStyle,'docked')
if ismac || ispc
matlab.graphics.internal.copyFigureHelper(fig_h);
else
%warning('Copy function to the clipboard only works if the figure is undocked.');
Helper.screencapture(fig_h,[],'clipboard');
end
else
pos = fig_h.Position;
Helper.screencapture(fig_h,[],'clipboard','position',[7,7,pos(3)-2,pos(4)]);
end
end

92
MOT Simulator/+Helper/ode5.m Normal file
View File

@ -0,0 +1,92 @@
function Y = ode5(odefun,tspan,y0,varargin)
%ODE5 Solve differential equations with a non-adaptive method of order 5.
% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates
% the system of differential equations y' = f(t,y) by stepping from T0 to
% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
% The vector Y0 is the initial conditions at T0. Each row in the solution
% array Y corresponds to a time specified in TSPAN.
%
% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters
% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).
%
% This is a non-adaptive solver. The step sequence is determined by TSPAN
% but the derivative function ODEFUN is evaluated multiple times per step.
% The solver implements the Dormand-Prince method of order 5 in a general
% framework of explicit Runge-Kutta methods.
%
% Example
% tspan = 0:0.1:20;
% y = ode5(@vdp1,tspan,[2 0]);
% plot(tspan,y(:,1));
% solves the system y' = vdp1(t,y) with a constant step size of 0.1,
% and plots the first component of the solution.
if ~isnumeric(tspan)
error('TSPAN should be a vector of integration steps.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
error(msg);
end
y0 = y0(:); % Make a column vector.
if ~isequal(size(y0),size(f0))
error('Inconsistent sizes of Y0 and f(t0,y0).');
end
neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
% Method coefficients -- Butcher's tableau
%
% C | A
% --+---
% | B
C = [1/5; 3/10; 4/5; 8/9; 1];
A = [ 1/5, 0, 0, 0, 0
3/40, 9/40, 0, 0, 0
44/45 -56/15, 32/9, 0, 0
19372/6561, -25360/2187, 64448/6561, -212/729, 0
9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];
B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
% More convenient storage
A = A.';
B = B(:);
nstages = length(B);
F = zeros(neq,nstages);
Y(:,1) = y0;
for i = 2:N
ti = tspan(i-1);
hi = h(i-1);
yi = Y(:,i-1);
% General explicit Runge-Kutta framework
F(:,1) = feval(odefun,ti,yi,varargin{:});
for stage = 2:nstages
tstage = ti + C(stage-1)*hi;
ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));
F(:,stage) = feval(odefun,tstage,ystage,varargin{:});
end
Y(:,i) = yi + F*(hi*B);
end
Y = Y.';

55
MOT Simulator/+Helper/onenoteccdata.m Normal file
View File

@ -0,0 +1,55 @@
cmap = zeros(16,16,3);
cmap(:,:,1) = [0.0000 0.0118 0.4510 0.0039 0.2078 0.1569 0.4078 0.4431 0.4510 0.1922 0.4235 0.4196 0.2235 0.4235 0.4039 0.4392
0.4471 0.1647 0.4157 0.0000 0.0235 0.4353 0.0314 0.4314 0.0196 0.2392 0.0667 0.0392 0.4431 0.3804 0.2941 0.4275
0.3686 0.3608 0.2000 0.2824 0.3059 0.0549 0.1804 0.1882 0.4392 0.4314 0.3255 0.0078 0.0902 0.1961 0.4353 0.1412
0.2314 0.3647 0.0353 0.3804 0.1647 0.2431 0.1686 0.2745 0.2980 0.4235 0.3922 0.4157 0.2784 0.3333 0.2510 0.0588
0.1020 0.0745 0.2549 0.0471 0.1216 0.4000 0.3961 0.2627 0.1098 0.1725 0.3098 0.4314 0.3529 0.3412 0.0784 0.0824
0.4471 0.1490 0.1804 0.3529 0.2196 0.3137 0.3255 0.0941 0.0078 0.3294 0.3765 0.2706 0.0510 0.0157 0.4275 0.1176
0.1294 0.1333 0.1725 0.3451 0.2118 0.3843 0.1255 0.1569 0.2118 0.1608 0.0353 0.2039 0.1608 0.4510 1.0000 0.8000
0.9882 0.6510 0.9961 0.4549 0.4549 0.6824 0.7882 0.5686 0.5373 0.5490 0.7765 0.7137 0.8510 0.7176 0.5020 0.4902
0.8941 0.9020 0.4745 0.8980 0.9098 0.4824 0.6471 0.6353 0.9922 0.9647 0.6353 0.4588 0.9647 0.9020 0.4980 0.8118
0.5059 0.4941 0.9686 0.4863 0.5451 0.9725 0.8980 0.5451 0.5333 0.6824 0.4588 0.8196 0.8314 0.8980 0.8941 0.9961
0.5255 0.8392 0.9804 0.5216 0.8588 0.8078 0.5176 0.7647 0.5608 0.9725 0.9059 0.4627 0.9882 0.8275 0.7725 0.8745
0.8235 0.8431 0.7373 1.0000 0.5137 0.4706 0.4784 0.7412 0.8863 0.9373 0.5529 0.5804 0.4510 0.9255 0.8235 0.8667
0.7569 0.8824 0.5294 0.5176 0.5373 0.9569 0.5294 0.4824 0.5098 0.5137 0.5569 0.8471 0.5098 0.9490 0.8706 0.9412
0.4902 0.6000 0.6980 0.7882 0.5490 0.7216 0.6431 0.4824 0.5569 0.4667 0.6627 0.9922 0.7804 0.8039 0.6275 0.7333
0.5725 0.5647 0.8549 0.7529 0.6235 0.8784 0.5922 0.7294 0.6118 0.7922 0.7843 0.6667 0.9294 0.6902 0.6784 0.9176
0.6706 0.7490 0.7961 0.5882 0.8627 0.4627 0.6196 0.7059 0.6078 0.9765 0.6549 0.6863 0.5373 0.7098 0.7176 0.7765];
cmap(:,:,2) = [0.0000 0.0078 0.2157 0.0000 0.0980 0.0745 0.1922 0.2157 0.2157 0.0902 0.2000 0.1961 0.1059 0.2039 0.1882 0.2078
0.2078 0.0784 0.2000 0.0000 0.0118 0.2118 0.0157 0.2039 0.0078 0.1137 0.0314 0.0196 0.2118 0.1804 0.1373 0.2078
0.1765 0.1725 0.0941 0.1333 0.1451 0.0275 0.0863 0.0902 0.2078 0.2078 0.1529 0.0039 0.0431 0.0941 0.2039 0.0667
0.1098 0.1725 0.0157 0.1804 0.0784 0.1137 0.0824 0.1333 0.1412 0.2000 0.1882 0.2000 0.1333 0.1569 0.1176 0.0275
0.0471 0.0353 0.1216 0.0196 0.0588 0.1922 0.1882 0.1255 0.0510 0.0824 0.1451 0.2039 0.1686 0.1647 0.0392 0.0392
0.2157 0.0706 0.0863 0.1686 0.1020 0.1490 0.1529 0.0431 0.0039 0.1569 0.1804 0.1255 0.0235 0.0078 0.2000 0.0549
0.0627 0.0627 0.0824 0.1647 0.1020 0.1843 0.0588 0.0745 0.1020 0.0784 0.0157 0.0980 0.0784 0.2157 1.0000 0.7137
0.9843 0.4980 0.9961 0.2235 0.2196 0.5412 0.6980 0.3843 0.3373 0.3569 0.6824 0.5922 0.7843 0.6000 0.2902 0.2706
0.8510 0.8588 0.2471 0.8549 0.8667 0.2627 0.4980 0.4784 0.9843 0.9490 0.4745 0.2235 0.9451 0.8627 0.2824 0.7333
0.2941 0.2784 0.9529 0.2667 0.3490 0.9569 0.8510 0.3490 0.3333 0.5451 0.2275 0.7412 0.7608 0.8549 0.8471 0.9922
0.3255 0.7686 0.9725 0.3176 0.8000 0.7255 0.3098 0.6627 0.3725 0.9647 0.8627 0.2314 0.9804 0.7529 0.6745 0.8235
0.7451 0.7765 0.6235 0.9961 0.3020 0.2431 0.2510 0.6314 0.8392 0.9098 0.3608 0.4000 0.2196 0.8902 0.7490 0.8078
0.6549 0.8353 0.3294 0.3137 0.3412 0.9373 0.3255 0.2588 0.2980 0.3059 0.3686 0.7843 0.3020 0.9255 0.8157 0.9176
0.2745 0.4275 0.5686 0.6980 0.3569 0.6039 0.4863 0.2627 0.3647 0.2392 0.5137 0.9922 0.6863 0.7216 0.4706 0.6196
0.3882 0.3765 0.7882 0.6471 0.4588 0.8275 0.4157 0.6118 0.4431 0.7059 0.6902 0.5255 0.8980 0.5569 0.5412 0.8824
0.5333 0.6392 0.7098 0.4078 0.8039 0.2314 0.4549 0.5804 0.4392 0.9647 0.5059 0.5529 0.3373 0.5882 0.5961 0.6784];
cmap(:,:,3) = [0.0000 0.0157 0.4980 0.0039 0.2314 0.1725 0.4627 0.5020 0.5020 0.2196 0.4745 0.4706 0.2510 0.4784 0.4510 0.4980
0.4941 0.1882 0.4667 0.0000 0.0275 0.4941 0.0353 0.4902 0.0196 0.2667 0.0745 0.0471 0.4902 0.4314 0.3294 0.4784
0.4196 0.4000 0.2235 0.3216 0.3412 0.0627 0.2039 0.2118 0.4863 0.4863 0.3608 0.0078 0.1020 0.2196 0.4824 0.1569
0.2588 0.4118 0.0392 0.4235 0.1843 0.2745 0.1882 0.3059 0.3373 0.4784 0.4392 0.4627 0.3137 0.3765 0.2824 0.0667
0.1137 0.0824 0.2863 0.0510 0.1373 0.4510 0.4471 0.2941 0.1216 0.1961 0.3490 0.4824 0.3961 0.3804 0.0902 0.0941
0.4980 0.1647 0.2000 0.4000 0.2431 0.3529 0.3647 0.1059 0.0118 0.3686 0.4196 0.3020 0.0549 0.0196 0.4824 0.1294
0.1451 0.1529 0.1922 0.3882 0.2392 0.4353 0.1412 0.1765 0.2353 0.1804 0.0353 0.2275 0.1843 0.5059 1.0000 0.8196
0.9882 0.6863 0.9961 0.5098 0.5098 0.7137 0.8118 0.6118 0.5843 0.5922 0.8000 0.7412 0.8627 0.7451 0.5529 0.5412
0.9059 0.9137 0.5255 0.9098 0.9176 0.5333 0.6824 0.6706 0.9922 0.9686 0.6706 0.5098 0.9647 0.9137 0.5490 0.8314
0.5569 0.5451 0.9725 0.5373 0.5922 0.9725 0.9059 0.5882 0.5804 0.7137 0.5137 0.8353 0.8510 0.9059 0.9020 0.9961
0.5725 0.8549 0.9843 0.5725 0.8745 0.8275 0.5647 0.7882 0.6039 0.9765 0.9137 0.5176 0.9882 0.8431 0.7961 0.8863
0.8392 0.8588 0.7647 1.0000 0.5608 0.5216 0.5294 0.7686 0.8980 0.9412 0.6000 0.6235 0.5059 0.9333 0.8431 0.8784
0.7804 0.8941 0.5765 0.5686 0.5843 0.9608 0.5765 0.5333 0.5569 0.5647 0.6039 0.8627 0.5608 0.9569 0.8863 0.9490
0.5412 0.6392 0.7294 0.8078 0.5961 0.7490 0.6784 0.5373 0.6000 0.5216 0.6941 0.9922 0.8039 0.8235 0.6667 0.7608
0.6157 0.6078 0.8667 0.7765 0.6588 0.8902 0.6314 0.7569 0.6510 0.8157 0.8039 0.7020 0.9373 0.7216 0.7098 0.9255
0.7059 0.7725 0.8196 0.6314 0.8784 0.5137 0.6549 0.7373 0.6471 0.9804 0.6902 0.7176 0.5804 0.7412 0.7451 0.8000];
%%
[cdata, cmap] = imread('onenote.png');

148
MOT Simulator/+Helper/parforNotifications.m Normal file
View File

@ -0,0 +1,148 @@
% Copyright (c) 2019 Andrea Alberti
%
% All rights reserved.
classdef parforNotifications < handle
properties
N; % number of iterations
text = 'Please wait ...'; % text to show
width = 50;
showWarning = true;
end
properties (GetAccess = public, SetAccess = private)
n;
end
properties (Access = private)
inProgress = false;
percent;
DataQueue;
usePercent;
Nstr;
NstrL;
lastComment;
end
methods
function this = parforNotifications()
this.DataQueue = parallel.pool.DataQueue;
afterEach(this.DataQueue, @this.updateStatus);
end
% Start progress bar
function PB_start(this,N,varargin)
assert(isscalar(N) && isnumeric(N) && N == floor(N) && N>0, 'Error: ''N'' must be a scalar positive integer.');
this.N = N;
p = inputParser;
addParameter(p,'message','Please wait: ');
addParameter(p,'usePercentage',true);
parse(p,varargin{:});
this.text = p.Results.message;
assert(ischar(this.text), 'Error: ''Message'' must be a string.');
this.usePercent = p.Results.usePercentage;
assert(isscalar(this.usePercent) && islogical(this.usePercent), 'Error: ''usePercentage'' must be a logical scalar.');
this.percent = 0;
this.n = 0;
this.lastComment = '';
if this.usePercent
fprintf('%s [%s]: %3d%%\n',this.text, char(32*ones(1,this.width)),0);
else
this.Nstr = sprintf('%d',this.N);
this.NstrL = numel(this.Nstr);
fprintf('%s [%s]: %s/%s\n',this.text, char(32*ones(1,this.width)),[char(32*ones(1,this.NstrL-1)),'0'],this.Nstr);
end
this.inProgress = true;
end
% Iterate progress bar
function PB_iterate(this,str)
if nargin == 1
send(this.DataQueue,'');
else
send(this.DataQueue,str);
end
end
function warning(this,warn_id,msg)
if this.showWarning
msg = struct('Action','Warning','Id',warn_id,'Message',msg);
send(this.DataQueue,msg);
end
end
function PB_reprint(this)
p = round(100*this.n/this.N);
this.percent = p;
cursor_pos=1+round((this.width-1)*p/100);
if p < 100
sep_char = '|';
else
sep_char = '.';
end
if this.usePercent
fprintf('%s [%s%s%s]: %3d%%\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p);
else
nstr=sprintf('%d',this.n);
fprintf('%s [%s%s%s]: %s/%s\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr);
end
end
function updateStatus(this,data)
if ischar(data)
this.n = this.n + 1;
p = round(100*this.n/this.N);
if p >= this.percent+1 || this.n == this.N
this.percent = p;
cursor_pos=1+round((this.width-1)*p/100);
if p < 100
sep_char = '|';
else
sep_char = '.';
end
if ~isempty(data)
comment = [' (',data,')'];
else
comment = '';
end
if this.usePercent
fprintf('%s%s%s%s]: %3d%%%s\n',char(8*ones(1,58+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p,comment);
else
nstr=sprintf('%d',this.n);
fprintf('%s%s%s%s]: %s/%s%s\n',char(8*ones(1,55+2*numel(this.Nstr)+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr,comment)
end
this.lastComment = comment;
if p == 100
this.inProgress = false;
end
end
else
switch data.Action
case 'Warning'
warning(data.Id,[data.Message,newline]);
if this.inProgress
this.PB_reprint();
end
end
end
end
end
end

820
MOT Simulator/+Helper/screencapture.m Normal file
View File

@ -0,0 +1,820 @@
function imageData = screencapture(varargin)
% screencapture - get a screen-capture of a figure frame, component handle, or screen area rectangle
%
% ScreenCapture gets a screen-capture of any Matlab GUI handle (including desktop,
% figure, axes, image or uicontrol), or a specified area rectangle located relative to
% the specified handle. Screen area capture is possible by specifying the root (desktop)
% handle (=0). The output can be either to an image file or to a Matlab matrix (useful
% for displaying via imshow() or for further processing) or to the system clipboard.
% This utility also enables adding a toolbar button for easy interactive screen-capture.
%
% Syntax:
% imageData = screencapture(handle, position, target, 'PropName',PropValue, ...)
%
% Input Parameters:
% handle - optional handle to be used for screen-capture origin.
% If empty/unsupplied then current figure (gcf) will be used.
% position - optional position array in pixels: [x,y,width,height].
% If empty/unsupplied then the handle's position vector will be used.
% If both handle and position are empty/unsupplied then the position
% will be retrieved via interactive mouse-selection.
% If handle is an image, then position is in data (not pixel) units, so the
% captured region remains the same after figure/axes resize (like imcrop)
% target - optional filename for storing the screen-capture, or the
% 'clipboard'/'printer' strings.
% If empty/unsupplied then no output to file will be done.
% The file format will be determined from the extension (JPG/PNG/...).
% Supported formats are those supported by the imwrite function.
% 'PropName',PropValue -
% optional list of property pairs (e.g., 'target','myImage.png','pos',[10,20,30,40],'handle',gca)
% PropNames may be abbreviated and are case-insensitive.
% PropNames may also be given in whichever order.
% Supported PropNames are:
% - 'handle' (default: gcf handle)
% - 'position' (default: gcf position array)
% - 'target' (default: '')
% - 'toolbar' (figure handle; default: gcf)
% this adds a screen-capture button to the figure's toolbar
% If this parameter is specified, then no screen-capture
% will take place and the returned imageData will be [].
%
% Output parameters:
% imageData - image data in a format acceptable by the imshow function
% If neither target nor imageData were specified, the user will be
% asked to interactively specify the output file.
%
% Examples:
% imageData = screencapture; % interactively select screen-capture rectangle
% imageData = screencapture(hListbox); % capture image of a uicontrol
% imageData = screencapture(0, [20,30,40,50]); % capture a small desktop region
% imageData = screencapture(gcf,[20,30,40,50]); % capture a small figure region
% imageData = screencapture(gca,[10,20,30,40]); % capture a small axes region
% imshow(imageData); % display the captured image in a matlab figure
% imwrite(imageData,'myImage.png'); % save the captured image to file
% img = imread('cameraman.tif');
% hImg = imshow(img);
% screencapture(hImg,[60,35,140,80]); % capture a region of an image
% screencapture(gcf,[],'myFigure.jpg'); % capture the entire figure into file
% screencapture(gcf,[],'clipboard'); % capture the entire figure into clipboard
% screencapture(gcf,[],'printer'); % print the entire figure
% screencapture('handle',gcf,'target','myFigure.jpg'); % same as previous, save to file
% screencapture('handle',gcf,'target','clipboard'); % same as previous, copy to clipboard
% screencapture('handle',gcf,'target','printer'); % same as previous, send to printer
% screencapture('toolbar',gcf); % adds a screen-capture button to gcf's toolbar
% screencapture('toolbar',[],'target','sc.bmp'); % same with default output filename
%
% Technical description:
% http://UndocumentedMatlab.com/blog/screencapture-utility/
%
% Bugs and suggestions:
% Please send to Yair Altman (altmany at gmail dot com)
%
% See also:
% imshow, imwrite, print
%
% Release history:
% 1.17 2016-05-16: Fix annoying warning about JavaFrame property becoming obsolete someday (yes, we know...)
% 1.16 2016-04-19: Fix for deployed application suggested by Dwight Bartholomew
% 1.10 2014-11-25: Added the 'print' target
% 1.9 2014-11-25: Fix for saving GIF files
% 1.8 2014-11-16: Fixes for R2014b
% 1.7 2014-04-28: Fixed bug when capturing interactive selection
% 1.6 2014-04-22: Only enable image formats when saving to an unspecified file via uiputfile
% 1.5 2013-04-18: Fixed bug in capture of non-square image; fixes for Win64
% 1.4 2013-01-27: Fixed capture of Desktop (root); enabled rbbox anywhere on desktop (not necesarily in a Matlab figure); enabled output to clipboard (based on Jiro Doke's imclipboard utility); edge-case fixes; added Java compatibility check
% 1.3 2012-07-23: Capture current object (uicontrol/axes/figure) if w=h=0 (e.g., by clicking a single point); extra input args sanity checks; fix for docked windows and image axes; include axes labels & ticks by default when capturing axes; use data-units position vector when capturing images; many edge-case fixes
% 1.2 2011-01-16: another performance boost (thanks to Jan Simon); some compatibility fixes for Matlab 6.5 (untested)
% 1.1 2009-06-03: Handle missing output format; performance boost (thanks to Urs); fix minor root-handle bug; added toolbar button option
% 1.0 2009-06-02: First version posted on <a href="http://www.mathworks.com/matlabcentral/fileexchange/authors/27420">MathWorks File Exchange</a>
% License to use and modify this code is granted freely to all interested, as long as the original author is
% referenced and attributed as such. The original author maintains the right to be solely associated with this work.
% Programmed and Copyright by Yair M. Altman: altmany(at)gmail.com
% $Revision: 1.17 $ $Date: 2016/05/16 17:59:36 $
% Ensure that java awt is enabled...
if ~usejava('awt')
error('YMA:screencapture:NeedAwt','ScreenCapture requires Java to run.');
end
% Ensure that our Java version supports the Robot class (requires JVM 1.3+)
try
robot = java.awt.Robot; %#ok<NASGU>
catch
uiwait(msgbox({['Your Matlab installation is so old that its Java engine (' version('-java') ...
') does not have a java.awt.Robot class. '], ' ', ...
'Without this class, taking a screen-capture is impossible.', ' ', ...
'So, either install JVM 1.3 or higher, or use a newer Matlab release.'}, ...
'ScreenCapture', 'warn'));
if nargout, imageData = []; end
return;
end
% Process optional arguments
paramsStruct = processArgs(varargin{:});
% If toolbar button requested, add it and exit
if ~isempty(paramsStruct.toolbar)
% Add the toolbar button
addToolbarButton(paramsStruct);
% Return the figure to its pre-undocked state (when relevant)
redockFigureIfRelevant(paramsStruct);
% Exit immediately (do NOT take a screen-capture)
if nargout, imageData = []; end
return;
end
% Convert position from handle-relative to desktop Java-based pixels
[paramsStruct, msgStr] = convertPos(paramsStruct);
% Capture the requested screen rectangle using java.awt.Robot
imgData = getScreenCaptureImageData(paramsStruct.position);
% Return the figure to its pre-undocked state (when relevant)
redockFigureIfRelevant(paramsStruct);
% Save image data in file or clipboard, if specified
if ~isempty(paramsStruct.target)
if strcmpi(paramsStruct.target,'clipboard')
if ~isempty(imgData)
imclipboard(imgData);
else
msgbox('No image area selected - not copying image to clipboard','ScreenCapture','warn');
end
elseif strncmpi(paramsStruct.target,'print',5) % 'print' or 'printer'
if ~isempty(imgData)
hNewFig = figure('visible','off');
imshow(imgData);
print(hNewFig);
delete(hNewFig);
else
msgbox('No image area selected - not printing screenshot','ScreenCapture','warn');
end
else % real filename
if ~isempty(imgData)
imwrite(imgData,paramsStruct.target);
else
msgbox(['No image area selected - not saving image file ' paramsStruct.target],'ScreenCapture','warn');
end
end
end
% Return image raster data to user, if requested
if nargout
imageData = imgData;
% If neither output formats was specified (neither target nor output data)
elseif isempty(paramsStruct.target) & ~isempty(imgData) %#ok ML6
% Ask the user to specify a file
%error('YMA:screencapture:noOutput','No output specified for ScreenCapture: specify the output filename and/or output data');
%format = '*.*';
formats = imformats;
for idx = 1 : numel(formats)
ext = sprintf('*.%s;',formats(idx).ext{:});
format(idx,1:2) = {ext(1:end-1), formats(idx).description}; %#ok<AGROW>
end
[filename,pathname] = uiputfile(format,'Save screen capture as');
if ~isequal(filename,0) & ~isequal(pathname,0) %#ok Matlab6 compatibility
try
filename = fullfile(pathname,filename);
imwrite(imgData,filename);
catch % possibly a GIF file that requires indexed colors
[imgData,map] = rgb2ind(imgData,256);
imwrite(imgData,map,filename);
end
else
% TODO - copy to clipboard
end
end
% Display msgStr, if relevant
if ~isempty(msgStr)
uiwait(msgbox(msgStr,'ScreenCapture'));
drawnow; pause(0.05); % time for the msgbox to disappear
end
return; % debug breakpoint
%% Process optional arguments
function paramsStruct = processArgs(varargin)
% Get the properties in either direct or P-V format
[regParams, pvPairs] = parseparams(varargin);
% Now process the optional P-V params
try
% Initialize
paramName = [];
paramsStruct = [];
paramsStruct.handle = [];
paramsStruct.position = [];
paramsStruct.target = '';
paramsStruct.toolbar = [];
paramsStruct.wasDocked = 0; % no false available in ML6
paramsStruct.wasInteractive = 0; % no false available in ML6
% Parse the regular (non-named) params in recption order
if ~isempty(regParams) & (isempty(regParams{1}) | ishandle(regParams{1}(1))) %#ok ML6
paramsStruct.handle = regParams{1};
regParams(1) = [];
end
if ~isempty(regParams) & isnumeric(regParams{1}) & (length(regParams{1}) == 4) %#ok ML6
paramsStruct.position = regParams{1};
regParams(1) = [];
end
if ~isempty(regParams) & ischar(regParams{1}) %#ok ML6
paramsStruct.target = regParams{1};
end
% Parse the optional param PV pairs
supportedArgs = {'handle','position','target','toolbar'};
while ~isempty(pvPairs)
% Disregard empty propNames (may be due to users mis-interpretting the syntax help)
while ~isempty(pvPairs) & isempty(pvPairs{1}) %#ok ML6
pvPairs(1) = [];
end
if isempty(pvPairs)
break;
end
% Ensure basic format is valid
paramName = '';
if ~ischar(pvPairs{1})
error('YMA:screencapture:invalidProperty','Invalid property passed to ScreenCapture');
elseif length(pvPairs) == 1
if isempty(paramsStruct.target)
paramsStruct.target = pvPairs{1};
break;
else
error('YMA:screencapture:noPropertyValue',['No value specified for property ''' pvPairs{1} '''']);
end
end
% Process parameter values
paramName = pvPairs{1};
if strcmpi(paramName,'filename') % backward compatibility
paramName = 'target';
end
paramValue = pvPairs{2};
pvPairs(1:2) = [];
idx = find(strncmpi(paramName,supportedArgs,length(paramName)));
if ~isempty(idx)
%paramsStruct.(lower(supportedArgs{idx(1)})) = paramValue; % incompatible with ML6
paramsStruct = setfield(paramsStruct, lower(supportedArgs{idx(1)}), paramValue); %#ok ML6
% If 'toolbar' param specified, then it cannot be left empty - use gcf
if strncmpi(paramName,'toolbar',length(paramName)) & isempty(paramsStruct.toolbar) %#ok ML6
paramsStruct.toolbar = getCurrentFig;
end
elseif isempty(paramsStruct.target)
paramsStruct.target = paramName;
pvPairs = {paramValue, pvPairs{:}}; %#ok (more readable this way, although a bit less efficient...)
else
supportedArgsStr = sprintf('''%s'',',supportedArgs{:});
error('YMA:screencapture:invalidProperty','%s \n%s', ...
'Invalid property passed to ScreenCapture', ...
['Supported property names are: ' supportedArgsStr(1:end-1)]);
end
end % loop pvPairs
catch
if ~isempty(paramName), paramName = [' ''' paramName '''']; end
error('YMA:screencapture:invalidProperty','Error setting ScreenCapture property %s:\n%s',paramName,lasterr); %#ok<LERR>
end
%end % processArgs
%% Convert position from handle-relative to desktop Java-based pixels
function [paramsStruct, msgStr] = convertPos(paramsStruct)
msgStr = '';
try
% Get the screen-size for later use
screenSize = get(0,'ScreenSize');
% Get the containing figure's handle
hParent = paramsStruct.handle;
if isempty(paramsStruct.handle)
paramsStruct.hFigure = getCurrentFig;
hParent = paramsStruct.hFigure;
else
paramsStruct.hFigure = ancestor(paramsStruct.handle,'figure');
end
% To get the acurate pixel position, the figure window must be undocked
try
if strcmpi(get(paramsStruct.hFigure,'WindowStyle'),'docked')
set(paramsStruct.hFigure,'WindowStyle','normal');
drawnow; pause(0.25);
paramsStruct.wasDocked = 1; % no true available in ML6
end
catch
% never mind - ignore...
end
% The figure (if specified) must be in focus
if ~isempty(paramsStruct.hFigure) & ishandle(paramsStruct.hFigure) %#ok ML6
isFigureValid = 1; % no true available in ML6
figure(paramsStruct.hFigure);
else
isFigureValid = 0; % no false available in ML6
end
% Flush all graphic events to ensure correct rendering
drawnow; pause(0.01);
% No handle specified
wasPositionGiven = 1; % no true available in ML6
if isempty(paramsStruct.handle)
% Set default handle, if not supplied
paramsStruct.handle = paramsStruct.hFigure;
% If position was not specified, get it interactively using RBBOX
if isempty(paramsStruct.position)
[paramsStruct.position, jFrameUsed, msgStr] = getInteractivePosition(paramsStruct.hFigure); %#ok<ASGLU> jFrameUsed is unused
paramsStruct.wasInteractive = 1; % no true available in ML6
wasPositionGiven = 0; % no false available in ML6
end
elseif ~ishandle(paramsStruct.handle)
% Handle was supplied - ensure it is a valid handle
error('YMA:screencapture:invalidHandle','Invalid handle passed to ScreenCapture');
elseif isempty(paramsStruct.position)
% Handle was supplied but position was not, so use the handle's position
paramsStruct.position = getPixelPos(paramsStruct.handle);
paramsStruct.position(1:2) = 0;
wasPositionGiven = 0; % no false available in ML6
elseif ~isnumeric(paramsStruct.position) | (length(paramsStruct.position) ~= 4) %#ok ML6
% Both handle & position were supplied - ensure a valid pixel position vector
error('YMA:screencapture:invalidPosition','Invalid position vector passed to ScreenCapture: \nMust be a [x,y,w,h] numeric pixel array');
end
% Capture current object (uicontrol/axes/figure) if w=h=0 (single-click in interactive mode)
if paramsStruct.position(3)<=0 | paramsStruct.position(4)<=0 %#ok ML6
%TODO - find a way to single-click another Matlab figure (the following does not work)
%paramsStruct.position = getPixelPos(ancestor(hittest,'figure'));
paramsStruct.position = getPixelPos(paramsStruct.handle);
paramsStruct.position(1:2) = 0;
paramsStruct.wasInteractive = 0; % no false available in ML6
wasPositionGiven = 0; % no false available in ML6
end
% First get the parent handle's desktop-based Matlab pixel position
parentPos = [0,0,0,0];
dX = 0;
dY = 0;
dW = 0;
dH = 0;
if ~isFigure(hParent)
% Get the reguested component's pixel position
parentPos = getPixelPos(hParent, 1); % no true available in ML6
% Axes position inaccuracy estimation
deltaX = 3;
deltaY = -1;
% Fix for images
if isImage(hParent) % | (isAxes(hParent) & strcmpi(get(hParent,'YDir'),'reverse')) %#ok ML6
% Compensate for resized image axes
hAxes = get(hParent,'Parent');
if all(get(hAxes,'DataAspectRatio')==1) % sanity check: this is the normal behavior
% Note 18/4/2013: the following fails for non-square images
%actualImgSize = min(parentPos(3:4));
%dX = (parentPos(3) - actualImgSize) / 2;
%dY = (parentPos(4) - actualImgSize) / 2;
%parentPos(3:4) = actualImgSize;
% The following should work for all types of images
actualImgSize = size(get(hParent,'CData'));
dX = (parentPos(3) - min(parentPos(3),actualImgSize(2))) / 2;
dY = (parentPos(4) - min(parentPos(4),actualImgSize(1))) / 2;
parentPos(3:4) = actualImgSize([2,1]);
%parentPos(3) = max(parentPos(3),actualImgSize(2));
%parentPos(4) = max(parentPos(4),actualImgSize(1));
end
% Fix user-specified img positions (but not auto-inferred ones)
if wasPositionGiven
% In images, use data units rather than pixel units
% Reverse the YDir
ymax = max(get(hParent,'YData'));
paramsStruct.position(2) = ymax - paramsStruct.position(2) - paramsStruct.position(4);
% Note: it would be best to use hgconvertunits, but:
% ^^^^ (1) it fails on Matlab 6, and (2) it doesn't accept Data units
%paramsStruct.position = hgconvertunits(hFig, paramsStruct.position, 'Data', 'pixel', hParent); % fails!
xLims = get(hParent,'XData');
yLims = get(hParent,'YData');
xPixelsPerData = parentPos(3) / (diff(xLims) + 1);
yPixelsPerData = parentPos(4) / (diff(yLims) + 1);
paramsStruct.position(1) = round((paramsStruct.position(1)-xLims(1)) * xPixelsPerData);
paramsStruct.position(2) = round((paramsStruct.position(2)-yLims(1)) * yPixelsPerData + 2*dY);
paramsStruct.position(3) = round( paramsStruct.position(3) * xPixelsPerData);
paramsStruct.position(4) = round( paramsStruct.position(4) * yPixelsPerData);
% Axes position inaccuracy estimation
if strcmpi(computer('arch'),'win64')
deltaX = 7;
deltaY = -7;
else
deltaX = 3;
deltaY = -3;
end
else % axes/image position was auto-infered (entire image)
% Axes position inaccuracy estimation
if strcmpi(computer('arch'),'win64')
deltaX = 6;
deltaY = -6;
else
deltaX = 2;
deltaY = -2;
end
dW = -2*dX;
dH = -2*dY;
end
end
%hFig = ancestor(hParent,'figure');
hParent = paramsStruct.hFigure;
elseif paramsStruct.wasInteractive % interactive figure rectangle
% Compensate for 1px rbbox inaccuracies
deltaX = 2;
deltaY = -2;
else % non-interactive figure
% Compensate 4px figure boundaries = difference betweeen OuterPosition and Position
deltaX = -1;
deltaY = 1;
end
%disp(paramsStruct.position) % for debugging
% Now get the pixel position relative to the monitor
figurePos = getPixelPos(hParent);
desktopPos = figurePos + parentPos;
% Now convert to Java-based pixels based on screen size
% Note: multiple monitors are automatically handled correctly, since all
% ^^^^ Java positions are relative to the main monitor's top-left corner
javaX = desktopPos(1) + paramsStruct.position(1) + deltaX + dX;
javaY = screenSize(4) - desktopPos(2) - paramsStruct.position(2) - paramsStruct.position(4) + deltaY + dY;
width = paramsStruct.position(3) + dW;
height = paramsStruct.position(4) + dH;
paramsStruct.position = round([javaX, javaY, width, height]);
%paramsStruct.position
% Ensure the figure is at the front so it can be screen-captured
if isFigureValid
figure(hParent);
drawnow;
pause(0.02);
end
catch
% Maybe root/desktop handle (root does not have a 'Position' prop so getPixelPos croaks
if isequal(double(hParent),0) % =root/desktop handle; handles case of hParent=[]
javaX = paramsStruct.position(1) - 1;
javaY = screenSize(4) - paramsStruct.position(2) - paramsStruct.position(4) - 1;
paramsStruct.position = [javaX, javaY, paramsStruct.position(3:4)];
end
end
%end % convertPos
%% Interactively get the requested capture rectangle
function [positionRect, jFrameUsed, msgStr] = getInteractivePosition(hFig)
msgStr = '';
try
% First try the invisible-figure approach, in order to
% enable rbbox outside any existing figure boundaries
f = figure('units','pixel','pos',[-100,-100,10,10],'HitTest','off');
drawnow; pause(0.01);
oldWarn = warning('off','MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame');
jf = get(handle(f),'JavaFrame');
warning(oldWarn);
try
jWindow = jf.fFigureClient.getWindow;
catch
try
jWindow = jf.fHG1Client.getWindow;
catch
jWindow = jf.getFigurePanelContainer.getParent.getTopLevelAncestor;
end
end
com.sun.awt.AWTUtilities.setWindowOpacity(jWindow,0.05); %=nearly transparent (not fully so that mouse clicks are captured)
jWindow.setMaximized(1); % no true available in ML6
jFrameUsed = 1; % no true available in ML6
msg = {'Mouse-click and drag a bounding rectangle for screen-capture ' ...
... %'or single-click any Matlab figure to capture the entire figure.' ...
};
catch
% Something failed, so revert to a simple rbbox on a visible figure
try delete(f); drawnow; catch, end %Cleanup...
jFrameUsed = 0; % no false available in ML6
msg = {'Mouse-click within any Matlab figure and then', ...
'drag a bounding rectangle for screen-capture,', ...
'or single-click to capture the entire figure'};
end
uiwait(msgbox(msg,'ScreenCapture'));
k = waitforbuttonpress; %#ok k is unused
%hFig = getCurrentFig;
%p1 = get(hFig,'CurrentPoint');
positionRect = rbbox;
%p2 = get(hFig,'CurrentPoint');
if jFrameUsed
jFrameOrigin = getPixelPos(f);
delete(f); drawnow;
try
figOrigin = getPixelPos(hFig);
catch % empty/invalid hFig handle
figOrigin = [0,0,0,0];
end
else
if isempty(hFig)
jFrameOrigin = getPixelPos(gcf);
else
jFrameOrigin = [0,0,0,0];
end
figOrigin = [0,0,0,0];
end
positionRect(1:2) = positionRect(1:2) + jFrameOrigin(1:2) - figOrigin(1:2);
if prod(positionRect(3:4)) > 0
msgStr = sprintf('%dx%d area captured',positionRect(3),positionRect(4));
end
%end % getInteractivePosition
%% Get current figure (even if its handle is hidden)
function hFig = getCurrentFig
oldState = get(0,'showHiddenHandles');
set(0,'showHiddenHandles','on');
hFig = get(0,'CurrentFigure');
set(0,'showHiddenHandles',oldState);
%end % getCurrentFig
%% Get ancestor figure - used for old Matlab versions that don't have a built-in ancestor()
function hObj = ancestor(hObj,type)
if ~isempty(hObj) & ishandle(hObj) %#ok for Matlab 6 compatibility
try
hObj = get(hObj,'Ancestor');
catch
% never mind...
end
try
%if ~isa(handle(hObj),type) % this is best but always returns 0 in Matlab 6!
%if ~isprop(hObj,'type') | ~strcmpi(get(hObj,'type'),type) % no isprop() in ML6!
try
objType = get(hObj,'type');
catch
objType = '';
end
if ~strcmpi(objType,type)
try
parent = get(handle(hObj),'parent');
catch
parent = hObj.getParent; % some objs have no 'Parent' prop, just this method...
end
if ~isempty(parent) % empty parent means root ancestor, so exit
hObj = ancestor(parent,type);
end
end
catch
% never mind...
end
end
%end % ancestor
%% Get position of an HG object in specified units
function pos = getPos(hObj,field,units)
% Matlab 6 did not have hgconvertunits so use the old way...
oldUnits = get(hObj,'units');
if strcmpi(oldUnits,units) % don't modify units unless we must!
pos = get(hObj,field);
else
set(hObj,'units',units);
pos = get(hObj,field);
set(hObj,'units',oldUnits);
end
%end % getPos
%% Get pixel position of an HG object - for Matlab 6 compatibility
function pos = getPixelPos(hObj,varargin)
persistent originalObj
try
stk = dbstack;
if ~strcmp(stk(2).name,'getPixelPos')
originalObj = hObj;
end
if isFigure(hObj) %| isAxes(hObj)
%try
pos = getPos(hObj,'OuterPosition','pixels');
else %catch
% getpixelposition is unvectorized unfortunately!
pos = getpixelposition(hObj,varargin{:});
% add the axes labels/ticks if relevant (plus a tiny margin to fix 2px label/title inconsistencies)
if isAxes(hObj) & ~isImage(originalObj) %#ok ML6
tightInsets = getPos(hObj,'TightInset','pixel');
pos = pos + tightInsets.*[-1,-1,1,1] + [-1,1,1+tightInsets(1:2)];
end
end
catch
try
% Matlab 6 did not have getpixelposition nor hgconvertunits so use the old way...
pos = getPos(hObj,'Position','pixels');
catch
% Maybe the handle does not have a 'Position' prop (e.g., text/line/plot) - use its parent
pos = getPixelPos(get(hObj,'parent'),varargin{:});
end
end
% Handle the case of missing/invalid/empty HG handle
if isempty(pos)
pos = [0,0,0,0];
end
%end % getPixelPos
%% Adds a ScreenCapture toolbar button
function addToolbarButton(paramsStruct)
% Ensure we have a valid toolbar handle
hFig = ancestor(paramsStruct.toolbar,'figure');
if isempty(hFig)
error('YMA:screencapture:badToolbar','the ''Toolbar'' parameter must contain a valid GUI handle');
end
set(hFig,'ToolBar','figure');
hToolbar = findall(hFig,'type','uitoolbar');
if isempty(hToolbar)
error('YMA:screencapture:noToolbar','the ''Toolbar'' parameter must contain a figure handle possessing a valid toolbar');
end
hToolbar = hToolbar(1); % just in case there are several toolbars... - use only the first
% Prepare the camera icon
icon = ['3333333333333333'; ...
'3333333333333333'; ...
'3333300000333333'; ...
'3333065556033333'; ...
'3000000000000033'; ...
'3022222222222033'; ...
'3022220002222033'; ...
'3022203110222033'; ...
'3022201110222033'; ...
'3022204440222033'; ...
'3022220002222033'; ...
'3022222222222033'; ...
'3000000000000033'; ...
'3333333333333333'; ...
'3333333333333333'; ...
'3333333333333333'];
cm = [ 0 0 0; ... % black
0 0.60 1; ... % light blue
0.53 0.53 0.53; ... % light gray
NaN NaN NaN; ... % transparent
0 0.73 0; ... % light green
0.27 0.27 0.27; ... % gray
0.13 0.13 0.13]; % dark gray
cdata = ind2rgb(uint8(icon-'0'),cm);
% If the button does not already exit
hButton = findall(hToolbar,'Tag','ScreenCaptureButton');
tooltip = 'Screen capture';
if ~isempty(paramsStruct.target)
tooltip = [tooltip ' to ' paramsStruct.target];
end
if isempty(hButton)
% Add the button with the icon to the figure's toolbar
hButton = uipushtool(hToolbar, 'CData',cdata, 'Tag','ScreenCaptureButton', 'TooltipString',tooltip, 'ClickedCallback',['screencapture(''' paramsStruct.target ''')']); %#ok unused
else
% Otherwise, simply update the existing button
set(hButton, 'CData',cdata, 'Tag','ScreenCaptureButton', 'TooltipString',tooltip, 'ClickedCallback',['screencapture(''' paramsStruct.target ''')']);
end
%end % addToolbarButton
%% Java-get the actual screen-capture image data
function imgData = getScreenCaptureImageData(positionRect)
if isempty(positionRect) | all(positionRect==0) | positionRect(3)<=0 | positionRect(4)<=0 %#ok ML6
imgData = [];
else
% Use java.awt.Robot to take a screen-capture of the specified screen area
rect = java.awt.Rectangle(positionRect(1), positionRect(2), positionRect(3), positionRect(4));
robot = java.awt.Robot;
jImage = robot.createScreenCapture(rect);
% Convert the resulting Java image to a Matlab image
% Adapted for a much-improved performance from:
% http://www.mathworks.com/support/solutions/data/1-2WPAYR.html
h = jImage.getHeight;
w = jImage.getWidth;
%imgData = zeros([h,w,3],'uint8');
%pixelsData = uint8(jImage.getData.getPixels(0,0,w,h,[]));
%for i = 1 : h
% base = (i-1)*w*3+1;
% imgData(i,1:w,:) = deal(reshape(pixelsData(base:(base+3*w-1)),3,w)');
%end
% Performance further improved based on feedback from Urs Schwartz:
%pixelsData = reshape(typecast(jImage.getData.getDataStorage,'uint32'),w,h).';
%imgData(:,:,3) = bitshift(bitand(pixelsData,256^1-1),-8*0);
%imgData(:,:,2) = bitshift(bitand(pixelsData,256^2-1),-8*1);
%imgData(:,:,1) = bitshift(bitand(pixelsData,256^3-1),-8*2);
% Performance even further improved based on feedback from Jan Simon:
pixelsData = reshape(typecast(jImage.getData.getDataStorage, 'uint8'), 4, w, h);
imgData = cat(3, ...
transpose(reshape(pixelsData(3, :, :), w, h)), ...
transpose(reshape(pixelsData(2, :, :), w, h)), ...
transpose(reshape(pixelsData(1, :, :), w, h)));
end
%end % getInteractivePosition
%% Return the figure to its pre-undocked state (when relevant)
function redockFigureIfRelevant(paramsStruct)
if paramsStruct.wasDocked
try
set(paramsStruct.hFigure,'WindowStyle','docked');
%drawnow;
catch
% never mind - ignore...
end
end
%end % redockFigureIfRelevant
%% Copy screen-capture to the system clipboard
% Adapted from http://www.mathworks.com/matlabcentral/fileexchange/28708-imclipboard/content/imclipboard.m
function imclipboard(imgData)
% Import necessary Java classes
import java.awt.Toolkit.*
import java.awt.image.BufferedImage
import java.awt.datatransfer.DataFlavor
% Add the necessary Java class (ImageSelection) to the Java classpath
if ~exist('ImageSelection', 'class')
% Obtain the directory of the executable (or of the M-file if not deployed)
%javaaddpath(fileparts(which(mfilename)), '-end');
if isdeployed % Stand-alone mode.
[status, result] = system('path'); %#ok<ASGLU>
MatLabFilePath = char(regexpi(result, 'Path=(.*?);', 'tokens', 'once'));
else % MATLAB mode.
MatLabFilePath = fileparts(mfilename('fullpath'));
end
javaaddpath(MatLabFilePath, '-end');
end
% Get System Clipboard object (java.awt.Toolkit)
cb = getDefaultToolkit.getSystemClipboard; % can't use () in ML6!
% Get image size
ht = size(imgData, 1);
wd = size(imgData, 2);
% Convert to Blue-Green-Red format
imgData = imgData(:, :, [3 2 1]);
% Convert to 3xWxH format
imgData = permute(imgData, [3, 2, 1]);
% Append Alpha data (not used)
imgData = cat(1, imgData, 255*ones(1, wd, ht, 'uint8'));
% Create image buffer
imBuffer = BufferedImage(wd, ht, BufferedImage.TYPE_INT_RGB);
imBuffer.setRGB(0, 0, wd, ht, typecast(imgData(:), 'int32'), 0, wd);
% Create ImageSelection object
% % custom java class
imSelection = ImageSelection(imBuffer);
% Set clipboard content to the image
cb.setContents(imSelection, []);
%end %imclipboard
%% Is the provided handle a figure?
function flag = isFigure(hObj)
flag = isa(handle(hObj),'figure') | isa(hObj,'matlab.ui.Figure');
%end %isFigure
%% Is the provided handle an axes?
function flag = isAxes(hObj)
flag = isa(handle(hObj),'axes') | isa(hObj,'matlab.graphics.axis.Axes');
%end %isFigure
%% Is the provided handle an image?
function flag = isImage(hObj)
flag = isa(handle(hObj),'image') | isa(hObj,'matlab.graphics.primitive.Image');
%end %isFigure
%%%%%%%%%%%%%%%%%%%%%%%%%% TODO %%%%%%%%%%%%%%%%%%%%%%%%%
% find a way in interactive-mode to single-click another Matlab figure for screen-capture

73
MOT Simulator/+Plotter/plotAngularDistributionForDifferentBeta.m Normal file
View File

@ -0,0 +1,73 @@
function plotAngularDistributionForDifferentBeta(obj, Beta)
f_h = Helper.getFigureByTag('AngDistForBeta');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Beta dependence';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
hold on
SimulationBeta = obj.Beta;
if ~ismember(SimulationBeta, Beta)
theta = linspace(0.0001,pi/2,1000);
J = zeros(1,length(theta));
for j=1:length(theta)
J(j) = obj.angularDistributionFunction(theta(j));
end
Norm = 0;
for j=1:length(J)
Norm = Norm+J(j)*sin(theta(j))*(theta(2)-theta(1));
end
J = J ./Norm*2;
J = J ./max(J);
plot([-flip(theta),theta], [flip(J),J],'DisplayName', ['\beta = ',num2str(SimulationBeta, '%.3f')], 'Linewidth',1.5)
end
for i=1:length(Beta)
theta = linspace(0.0001,pi/2,1000);
J = zeros(1,length(theta));
obj.NozzleLength = (2 * obj.NozzleRadius) / Beta(i);
for j=1:length(theta)
J(j) = obj.angularDistributionFunction(theta(j));
end
Norm = 0;
for j=1:length(J)
Norm = Norm+J(j)*sin(theta(j))*(theta(2)-theta(1));
end
J = J ./Norm*2;
J = J ./max(J);
if Beta(i) ~= SimulationBeta
plot([-flip(theta),theta], [flip(J),J],'DisplayName',['\beta = ',num2str(Beta(i))], 'LineStyle', '--', 'Linewidth',1.5)
else
plot([-flip(theta),theta], [flip(J),J],'DisplayName',['\beta = ',num2str(Beta(i))], 'Linewidth',1.5)
end
end
hold off
leg = legend;
hXLabel = xlabel('\theta (rad)');
hYLabel = ylabel('J(\theta)');
hTitle = sgtitle('Angular Distribution (Transition Flow)');
set([hXLabel, hYLabel, leg] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
grid on
Helper.bringFiguresWithTagInForeground();
end

37
MOT Simulator/+Plotter/plotCaptureVelocityVsAngle.m Normal file
View File

@ -0,0 +1,37 @@
function plotCaptureVelocityVsAngle(OvenObj, MOTObj)
theta = linspace(0, OvenObj.ExitDivergence, 1000);
CaptureVelocity = zeros(length(theta),3);
for i=1:length(theta)
CaptureVelocity(i,:) = MOTObj.calculateCaptureVelocity(OvenObj, [-OvenObj.OvenDistance,0,0], [cos(theta(i)),0,sin(theta(i))]);
end
f_h = Helper.getFigureByTag('Capture Velocity');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Capture Velocity';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
plot(theta .* 1e+03, sqrt(CaptureVelocity(:,1).^2+CaptureVelocity(:,2).^2+CaptureVelocity(:,3).^2), 'Linewidth', 1.5)
hXLabel = xlabel('\theta (mrad)');
hYLabel = ylabel('Velocity (m/s)');
hTitle = sgtitle('Capture velocity for different angles of divergence');
set([hXLabel, hYLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 14 );
grid on
Helper.bringFiguresWithTagInForeground();
end

18
MOT Simulator/+Plotter/plotConfidenceIntervalRegion.m Normal file
View File

@ -0,0 +1,18 @@
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

85
MOT Simulator/+Plotter/plotDynamicalQuantities.m Normal file
View File

@ -0,0 +1,85 @@
function plotDynamicalQuantities(OvenObj, MOTObj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, varargin)
p = inputParser;
p.KeepUnmatched = true;
addRequired(p, 'OvenObject' , @isobject);
addRequired(p, 'MOTObject' , @isobject);
addRequired(p, 'MaximumVelocity' , @(x) assert(isnumeric(x) && isscalar(x)));
addRequired(p, 'IncidentAtomDirection' , @(x) assert(isnumeric(x) && isscalar(x)));
addRequired(p, 'IncidentAtomPosition' , @(x) assert(isnumeric(x) && isscalar(x)));
addParameter(p, 'PlotPositions' , false, @islogical);
addParameter(p, 'PlotVelocities' , false, @islogical);
p.parse(OvenObj, MOTObj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, varargin{:})
MaximumVelocity = p.Results.MaximumVelocity;
IncidentAtomDirection = p.Results.IncidentAtomDirection;
IncidentAtomPosition = p.Results.IncidentAtomDirection;
PlotPositions = p.Results.PlotPositions;
PlotVelocities = p.Results.PlotVelocities;
f_h = Helper.getFigureByTag('Trajectories Plot');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Trajectories Plot';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/50 screensize(4)/3.5] 1870 600];
N = MOTObj.NumberOfAtoms;
Theta = IncidentAtomDirection;
z = IncidentAtomPosition;
L = OvenObj.OvenDistance * 2;
T = MOTObj.SimulationTime;
tau = MOTObj.TimeStep;
Y = linspace(0,MaximumVelocity,N);
DynamicalQuantities = zeros(length(Y),int64(T/tau),6);
for i=1:length(Y)
x =-L/2;
vx = Y(i)*cos(Theta);
vz = Y(i)*sin(Theta);
r = [x,0,z];
v = [vx,0,vz];
DynamicalQuantities(i,:,:) = MOTObj.solver(r, v);
end
LabelStringArrayForPositions = {'x (mm)', 'y (mm)', 'z (mm)'};
LabelStringArrayForVelocities = {'v_x (m/s)', 'v_y (m/s)', 'v_z (m/s)'};
colors = { [0, 0.4470, 0.7410], [0.8500, 0.3250, 0.0980], [0.4940, 0.1840, 0.5560]};
t = linspace(0, T, T/tau) * 1e+3;
for i=1:length(Y)
for j = 1:3
if PlotPositions
subplot(1, 3, j)
hold on
plot(t, DynamicalQuantities(i,:,j) * 1e+3,'Color', colors{j},'linewidth',1.3)
hXLabel = xlabel('Time (ms)');
hYLabel = ylabel(LabelStringArrayForPositions{j});
hold off
set([hXLabel, hYLabel], 'FontSize', 14);
elseif PlotVelocities
subplot(1, 3, j)
hold on
plot(t, DynamicalQuantities(i,:,j+3),'Color', colors{j},'linewidth',1.3)
hXLabel = xlabel('Time (ms)');
hYLabel = ylabel(LabelStringArrayForVelocities{j});
hold off
set([hXLabel, hYLabel], 'FontSize', 14);
end
end
end
hold off
hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", MOTObj.MagneticGradient));
set(hTitle, 'FontSize', 18);
Helper.bringFiguresWithTagInForeground();
end

54
MOT Simulator/+Plotter/plotFreeMolecularFluxVsTemp.m Normal file
View File

@ -0,0 +1,54 @@
function plotFreeMolecularFluxVsTemp(obj, Temperature)
f_h = Helper.getFigureByTag('Tube Flux');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Tube Flux';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/4.5 screensize(4)/4.5] 920 750];
hold on
for i=1:length(Temperature)
beta = linspace(0.01,0.5,200);
obj.OvenTemperature = Temperature(i);
flux = zeros(1,length(beta));
for j=1:length(beta)
obj.NozzleLength = (2 * obj.NozzleRadius) / beta(j);
[ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor();
flux(j) = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux();
end
plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5)
end
set(gca,'yscale','log')
obj.restoreDefaults();
xline(obj.Beta, 'k--','Linewidth', 0.5);
fmf = obj.ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux();
yline(fmf, 'k--', 'Linewidth', 1.5);
textstring = [sprintf('%1.e',fmf) ' atoms/s for ' '$$ \beta $$ = ' num2str(obj.Beta, '%.2f') sprintf(' @ %.2f K', obj.OvenTemperatureinKelvin)];
txt = text((obj.Beta + 0.05*obj.Beta), (max(fmf) + 0.2*fmf), textstring, 'Interpreter','latex');
hold off
leg = legend('Location', 'southeast');
leg.String = leg.String(1:end-2); % Remove xline and yline legend entries
hXLabel = xlabel('\beta');
hYLabel = ylabel('Flux (atoms/s)');
hTitle = sgtitle('Total Flux from a Tube (Free Molecular Flow)');
set([hXLabel, hYLabel, txt, leg] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
grid on
Helper.bringFiguresWithTagInForeground();
end

72
MOT Simulator/+Plotter/plotInitialVeloctiySamplingVsAngle.m Normal file
View File

@ -0,0 +1,72 @@
function plotInitialVeloctiySamplingVsAngle(obj, MOTObj, NumberOfBins)
n = obj.NumberOfAtoms;
VelocitySamples = obj.initialVelocitySampling(MOTObj);
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.ExitDivergence,1000);
AngleBinSize = obj.ExitDivergence / 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('|v| (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.ExitDivergence.* 1e+03, 'k--', 'Linewidth', 1.5);
text((obj.ExitDivergence - 0.5 * obj.ExitDivergence).* 1e+03, max(h_1.Values) - 0.45 * max(h_1.Values), 'Maximum Allowed Divergence ----------->');
hXLabel_2 = xlabel('\theta (mrad)');
hYLabel_2 = ylabel('Counts');
hold off
hTitle = sgtitle('Sampled Velocities');
set([hXLabel_1, hXLabel_2, hYLabel_1, hYLabel_2] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

43
MOT Simulator/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m Normal file
View File

@ -0,0 +1,43 @@
function plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
TemperatureinKelvin = TemperatureinCelsius + 273.15;
Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+TemperatureinKelvin)); % Vapor Pressure Dysprosium for the given oven temperature
MeanFreepath = (Helper.PhysicsConstants.BoltzmannConstant*TemperatureinKelvin)./(sqrt(2) * ( pi * (2*281e-12)^2) * Dy164VapourPressure); %free path at above tempeture
f_h = Helper.getFigureByTag('Dysprosium Gas Properties');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Dysprosium Gas Properties';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
yyaxis left
semilogy(TemperatureinCelsius, Dy164VapourPressure, 'Color', '#0071BB', 'Linewidth',1.5);
hYLabel_1 = ylabel('Vapor Pressure (mbar)');
yyaxis right
semilogy(TemperatureinCelsius, MeanFreepath.*1000, 'Color', '#36B449', 'Linewidth',1.5)
hYLabel_2 = ylabel('Mean Free Path (mm)');
hXLabel = xlabel('Temperature (℃)');
ax = gca;
ax.YAxis(1).Color = '#0071BB';
ax.YAxis(2).Color = '#36B449';
hTitle = sgtitle('^{164}Dy Gas');
set([hXLabel, hYLabel_1, hYLabel_2] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
grid on
Helper.bringFiguresWithTagInForeground();
end

84
MOT Simulator/+Plotter/plotPhaseSpaceWithAccelerationField.m Normal file
View File

@ -0,0 +1,84 @@
function plotPhaseSpaceWithAccelerationField(OvenObj, MOTObj, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
f_h = Helper.getFigureByTag('Phase Space Plot');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Phase Space Plot';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
N = MOTObj.NumberOfAtoms;
L = OvenObj.OvenDistance * 2;
Theta = IncidentAtomDirection;
z = IncidentAtomPosition;
T = MOTObj.SimulationTime;
tau = MOTObj.TimeStep;
[X,Y] = meshgrid(-L/2:L/NumberOfBins:L/2,-MaximumVelocity:2*MaximumVelocity/NumberOfBins:MaximumVelocity);
a=zeros(NumberOfBins+1,NumberOfBins+1,3);
% MOTObj.restoreDefaults();
for i=1:length(X)
for j=1:length(Y)
a(i,j,:) = MOTObj.calculateTotalAcceleration([X(1,i), 0, z], [Y(j,1)*cos(Theta),0,Y(j,1)*sin(Theta)]);
end
end
for i=1:length(X)
for j=1:length(Y)
if isnan(a(i,j,1)) || isnan(a(i,j,2)) || isnan(a(i,j,3))
a(i,j,1)=0;
a(i,j,2)=0;
a(i,j,3)=0;
end
end
end
pcolor(X',Y',a(:,:,1))
hold on
col=colorbar;
col.Label.String='Aceleration (m/s^2)';
shading flat
%-------------------------------------------------------------------------
Y = linspace(MinimumVelocity, MaximumVelocity,N);
DynamicalQuantities = zeros(length(Y),int64(T/tau),6);
for i=1:length(Y)
x =-L/2;
vx = Y(i)*cos(Theta);
vz = Y(i)*sin(Theta);
r = [x,0,z];
v = [vx,0,vz];
DynamicalQuantities(i,:,:) = MOTObj.solver(r, v);
end
hold on
count = 0;
for i=1:length(Y)
if DynamicalQuantities(i,end,2) > 0
plot(DynamicalQuantities(i,:,1),DynamicalQuantities(i,:,4),'w','linewidth',1.3)
count = count + 1;
end
end
hold off
hXLabel = xlabel('Position: Along the x-axis (m)');
hYLabel = ylabel('Velocity (m/s)');
hTitle = sgtitle(sprintf("Magnetic gradient = %.3f T/m", MOTObj.MagneticGradient));
set([hXLabel, hYLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

55
MOT Simulator/+Plotter/plotPositionAndVelocitySampling.m Normal file
View File

@ -0,0 +1,55 @@
function plotPositionAndVelocitySampling(NumberOfBins, initialPositions, initialVelocities)
f_h = Helper.getFigureByTag('RejectionSampling');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Sampling';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/7 screensize(4)/7] 1.357e+03 770];
subplot(3,2,1)
histogram(initialPositions(:, 1)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component')
xlabel('Positions (mm)','FontSize', 14)
ylabel('Counts','FontSize', 14)
legend('FontSize', 14)
subplot(3,2,3)
histogram(initialPositions(:, 2)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component')
xlabel('Positions (mm)','FontSize', 14)
ylabel('Counts','FontSize', 14)
legend('FontSize', 14)
subplot(3,2,5)
histogram(initialPositions(:, 3)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component')
xlabel('Positions (mm)','FontSize', 14)
ylabel('Counts','FontSize', 14)
legend('FontSize', 14)
subplot(3,2,2)
histogram(initialVelocities(:, 1),NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component')
xlabel('Velocities (m/s)','FontSize', 14)
ylabel('Counts','FontSize', 14)
legend('FontSize', 14)
subplot(3,2,4)
histogram(initialVelocities(:, 2),NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component')
xlabel('Velocities (m/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)
ylabel('Counts','FontSize', 14)
legend('FontSize', 14)
sgtitle('Rejection sampling for initial distributions','FontSize', 18)
Helper.bringFiguresWithTagInForeground();
end

90
MOT Simulator/+Plotter/plotResultForOneParameterScan.m Normal file
View File

@ -0,0 +1,90 @@
function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
p = inputParser;
p.KeepUnmatched = true;
addRequired(p, 'ParameterArray', @isvector)
addRequired(p, 'QuantityOfInterestArray', @ismatrix)
addParameter(p, 'RescalingFactorForParameter', 1, @isscalar)
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)
p.parse(XParameter, YQuantity, varargin{:})
XParameter = p.Results.ParameterArray;
RescalingFactorForXParameter = p.Results.RescalingFactorForParameter;
XLabelString = p.Results.XLabelString;
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;
f_h = Helper.getFigureByTag('One-Parameter Scan');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'One-Parameter Scan';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
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')
else
plot(RescaledXParameter, RescaledYQuantity, '--o', 'Linewidth', 1.5);
end
if CIForYQuantity
RescaledCIArray = CIArray .* RescalingFactorForYQuantity;
Plotter.plotConfidenceIntervalRegion(RescaledXParameter, RescaledCIArray(:,1), RescaledCIArray(:,2));
end
hold off
xlim([0 inf])
ylim([0 inf])
hXLabel = xlabel(XLabelString);
hYLabel = ylabel(YLabelString);
hTitle = sgtitle(TitleString);
set([hXLabel, hYLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
grid on
Helper.bringFiguresWithTagInForeground();
end

82
MOT Simulator/+Plotter/plotResultForThreeParameterScan.m Normal file
View File

@ -0,0 +1,82 @@
function plotResultForThreeParameterScan(XParameter, YParameter, DeltaParameter, ZQuantity, varargin)
p = inputParser;
p.KeepUnmatched = true;
addRequired(p, 'FirstParameterArray', @isvector)
addRequired(p, 'SecondParameterArray', @isvector)
addRequired(p, 'ThirdParameterArray', @isvector)
addRequired(p, 'QuantityOfInterestArray', @ismatrix)
addParameter(p, 'RescalingFactorForFirstParameter', 1, @isscalar)
addParameter(p, 'XLabelString', 'X parameter', @ischar)
addParameter(p, 'RescalingFactorForSecondParameter', 1, @isscalar)
addParameter(p, 'YLabelString', 'Y parameter', @ischar)
addParameter(p, 'RescalingFactorForThirdParameter', 1, @isscalar)
addParameter(p, 'RescalingFactorForQuantityOfInterest', 1, @isscalar)
addParameter(p, 'ZLabelString', 'Z parameter', @ischar)
addParameter(p, 'PlotTitleString', '<third parameter value>', @ischar)
addParameter(p, 'FigureTitleString', 'Third-Parameter Scan', @ischar)
p.parse(XParameter, YParameter, DeltaParameter, ZQuantity, varargin{:})
XParameter = p.Results.FirstParameterArray;
RescalingFactorForXParameter = p.Results.RescalingFactorForFirstParameter;
XLabelString = p.Results.XLabelString;
YParameter = p.Results.SecondParameterArray;
RescalingFactorForYParameter = p.Results.RescalingFactorForSecondParameter;
YLabelString = p.Results.YLabelString;
DeltaParameter = p.Results.ThirdParameterArray;
RescalingFactorForDeltaParameter = p.Results.RescalingFactorForThirdParameter;
ZQuantity = p.Results.QuantityOfInterestArray;
RescalingFactorForZQuantity = p.Results.RescalingFactorForQuantityOfInterest;
ZLabelString = p.Results.ZLabelString;
PlotTitleString = p.Results.PlotTitleString;
FigureTitleString = p.Results.FigureTitleString;
f_h = Helper.getFigureByTag('Three-Parameter Scan');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Three-Parameter Scan';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/60 screensize(4)/10] 1530 870];
RescaledXParameter = XParameter .* RescalingFactorForXParameter;
RescaledYParameter = YParameter .* RescalingFactorForYParameter;
RescaledDeltaParameter = DeltaParameter .* RescalingFactorForDeltaParameter;
tiledlayout(ceil(length(DeltaParameter)/3), 3)
for i = 1:length(DeltaParameter)
nexttile
RescaledZQuantity = ZQuantity{i} .* RescalingFactorForZQuantity;
imagesc(RescaledXParameter, RescaledYParameter, RescaledZQuantity(:,:)');
set(gca,'YDir','normal');
hXLabel = xlabel(XLabelString);
hYLabel = ylabel(YLabelString);
hPlotLabel = title(sprintf(PlotTitleString, RescaledDeltaParameter(i)));
set([hXLabel, hYLabel, hPlotLabel] , ...
'FontSize' , 14);
end
caxis([min(min(min(RescaledZQuantity))) max(max(max(RescaledZQuantity)))]);
shading flat;
c = colorbar;
c.Label.String= ZLabelString;
c.Label.FontSize = 14;
c.Layout.Tile = 'east';
hTitle = sgtitle(FigureTitleString);
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

74
MOT Simulator/+Plotter/plotResultForTwoParameterScan.m Normal file
View File

@ -0,0 +1,74 @@
function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, varargin)
p = inputParser;
p.KeepUnmatched = true;
addRequired(p, 'FirstParameterArray', @isvector)
addRequired(p, 'SecondParameterArray', @isvector)
addRequired(p, 'QuantityOfInterestArray', @ismatrix)
addParameter(p, 'RescalingFactorForFirstParameter', 1, @isscalar)
addParameter(p, 'XLabelString', 'X parameter', @ischar)
addParameter(p, 'RescalingFactorForSecondParameter', 1, @isscalar)
addParameter(p, 'YLabelString', 'Y parameter', @ischar)
addParameter(p, 'RescalingFactorForQuantityOfInterest', 1, @isscalar)
addParameter(p, 'ZLabelString', 'Z parameter', @ischar)
addParameter(p, 'TitleString', 'Two-Parameter Scan', @ischar)
p.parse(XParameter, YParameter, ZQuantity, varargin{:})
XParameter = p.Results.FirstParameterArray;
RescalingFactorForXParameter = p.Results.RescalingFactorForFirstParameter;
XLabelString = p.Results.XLabelString;
YParameter = p.Results.SecondParameterArray;
RescalingFactorForYParameter = p.Results.RescalingFactorForSecondParameter;
YLabelString = p.Results.YLabelString;
ZQuantity = p.Results.QuantityOfInterestArray;
RescalingFactorForZQuantity = p.Results.RescalingFactorForQuantityOfInterest;
ZLabelString = p.Results.ZLabelString;
TitleString = p.Results.TitleString;
f_h = Helper.getFigureByTag('Two-Parameter Scan');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Two-Parameter Scan';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
RescaledXParameter = XParameter .* RescalingFactorForXParameter;
RescaledYParameter = YParameter .* RescalingFactorForYParameter;
RescaledZQuantity = ZQuantity .* RescalingFactorForZQuantity;
imagesc(RescaledXParameter, RescaledYParameter, RescaledZQuantity(:,:)');
% hold on
%
% contour(RescaledXParameter, RescaledYParameter, RescaledZQuantity(:,:)', 'Color', 'r', 'Linewidth', 4, 'ShowText','on')
set(gca,'YDir','normal');
caxis([min(min(min(RescaledZQuantity))) max(max(max(RescaledZQuantity)))]);
hXLabel = xlabel(XLabelString);
hYLabel = ylabel(YLabelString);
shading flat;
c = colorbar;
c.Label.String= ZLabelString;
c.Label.FontSize = 14;
hTitle = sgtitle(TitleString);
set([hXLabel, hYLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

72
MOT Simulator/+Plotter/visualizeMagneticField.m Normal file
View File

@ -0,0 +1,72 @@
function visualizeMagneticField(obj, x_range, y_range, z_range)
f_h = Helper.getFigureByTag('VisualizeMagneticFieldFor2DMOT');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Visualization';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 820 645];
xmin = x_range(1);
xmax = x_range(2);
ymin = y_range(1);
ymax = y_range(2);
zmin = z_range(1);
zmax = z_range(2);
dx = (xmax-xmin)/8;
dy = (ymax-ymin)/8;
dz = (zmax-zmin)/8;
if dx ~= 0
xm = xmin:dx:xmax;
else
xm = zeros(1,5);
end
if dy ~= 0
ym = ymin:dy:ymax;
else
ym = zeros(1,5);
end
if dz ~= 0
zm = zmin:dz:zmax;
else
zm = zeros(1,5);
end
[meshx,meshy,meshz] = meshgrid(xm,ym,zm); % construct data points
switch obj.SimulationMode
case '2D'
alpha = obj.MagneticGradient;
Bx = @(x,y,z) alpha .* z;
By = @(x,y,z) 0 .* y;
Bz = @(x,y,z) alpha .* x;
Bx_val = Bx(meshx, meshy, meshz);
By_val = By(meshx, meshy, meshz);
Bz_val = Bz(meshx, meshy, meshz);
case '3D'
% Development in progress
end
quiver3(meshx, meshy, meshz, Bx_val, By_val, Bz_val, 'Color', ' #6600ff');
axis equal
hXLabel = xlabel('x');
hYLabel = ylabel('y');
hZLabel = zlabel('z');
hTitle = sgtitle('Magnetic Field for 2-D MOT');
set([hXLabel, hYLabel, hZLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
end

236
MOT Simulator/+Scripts/optimizingForSidebandEnhancement.m Normal file
View File

@ -0,0 +1,236 @@
OptionsStruct = struct;
OptionsStruct.ErrorEstimationMethod = 'bootstrap'; % 'jackknife' | 'bootstrap'
OptionsStruct.NumberOfAtoms = 5000;
OptionsStruct.TimeStep = 50e-06; % in s
OptionsStruct.SimulationTime = 5e-03; % in s
OptionsStruct.SpontaneousEmission = true;
OptionsStruct.SidebandBeam = true;
OptionsStruct.PushBeam = true;
OptionsStruct.Gravity = true;
OptionsStruct.BackgroundCollision = true;
OptionsStruct.SaveData = true;
% OptionsStruct.SaveDirectory = '';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
Oven = Simulator.Oven(options{:});
MOT2D = Simulator.TwoDimensionalMOT(options{:});
Beams = MOT2D.Beams;
%%
MOT2D.NumberOfAtoms = 10000;
MOT2D.TotalPower = 0.8;
MOT2D.MagneticGradient = 0.4; 0;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Waist = 15e-03;
CoolingBeam.Detuning = -1.67*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Waist = 15e-03;
NumberOfPointsForFirstParam = 20; %iterations of the simulation
NumberOfPointsForSecondParam = 20;
DetuningArray = linspace(-1.0, -6.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth;
PowerArray = linspace(0, 0.8, NumberOfPointsForSecondParam) * MOT2D.TotalPower;
tStart = tic;
[LoadingRateArray, ~, ~] = Scripts.scanForSidebandEnhancement(Oven, MOT2D, 'Blue', 'BlueSideband', DetuningArray, PowerArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
if MOT2D.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
MOT2D.Results = LoadingRate;
SaveFolder = [MOT2D.SaveDirectory filesep 'Results'];
Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = MOT2D;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
MOT2D.SidebandBeam = false;
MOT2D.PushBeam = false;
CoolingBeam.Power = MOT2D.TotalPower;
[LoadingRate, ~] = MOT2D.runSimulation(Oven);
EnhancementFactorArray = LoadingRateArray ./ LoadingRate;
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString = 'Sideband Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Sideband Beam Power (mW)';
OptionsStruct.RescalingFactorForQuantityOfInterest = 1;
OptionsStruct.ZLabelString = 'Enhancement Factor (\eta)';
% OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForTwoParameterScan(DetuningArray, PowerArray, EnhancementFactorArray, options{:})
%% Magnetic gradient scan
MOT2D.NumberOfAtoms = 10000;
MOT2D.TotalPower = 0.4;
MOT2D.SidebandBeam = true;
NumberOfPointsForFirstParam = 10; %iterations of the simulation
NumberOfPointsForSecondParam = 10;
NumberOfPointsForThirdParam = 6;
DetuningArray = linspace(-0.5, -5.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth;
PowerArray = linspace(0.1, 1.0, NumberOfPointsForSecondParam) * MOT2D.TotalPower;
MagneticGradientArray = linspace(30, 50, NumberOfPointsForThirdParam) * 1e-02;
Beams = MOT2D.Beams;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
LoadingRateArray = {};
tStart = tic;
for k=1:NumberOfPointsForThirdParam
eval(sprintf('MOT2D.%s = %d;', 'MagneticGradient', MagneticGradientArray(k)));
lrmatrix = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
for i=1:NumberOfPointsForFirstParam
eval(sprintf('SidebandBeam.%s = %d;', 'Detuning', DetuningArray(i)));
for j=1:NumberOfPointsForSecondParam
eval(sprintf('SidebandBeam.%s = %d;', 'Power', PowerArray(j)));
eval(sprintf('CoolingBeam.%s = %d;', 'Power', MOT2D.TotalPower - PowerArray(j)));
[lrmatrix(i,j), ~, ~] = MOT2D.runSimulation(Oven);
end
end
LoadingRateArray{end+1} = lrmatrix;
end
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString = 'Sideband Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Sideband Beam Waist (mm)';
OptionsStruct.RescalingFactorForThirdParameter = 100;
OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-9;
OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)';
OptionsStruct.PlotTitleString = 'Magnetic Gradient = %.0f (G/cm)';
OptionsStruct.FigureTitleString = sprintf('Oven-2DMOT Distance = %.1f (mm); Total Beam Power = %d (mW)', Oven.OvenDistance * 1000, MOT2D.TotalPower*1000);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForThreeParameterScan(DetuningArray, PowerArray, MagneticGradientArray, LoadingRateArray, options{:})
clear OptionsStruct
%%
MOT2D.NumberOfAtoms = 10000;
MOT2D.TotalPower = 0.4;
MOT2D.MagneticGradient = 0.4; 0;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power = 0.2;
CoolingBeam.Detuning = -1.3*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power = 0.2;
NumberOfPointsForFirstParam = 20; %iterations of the simulation
NumberOfPointsForSecondParam = 20;
DetuningArray = linspace(-1.0, -6.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth;
BeamWaistArray = linspace(10, 25, NumberOfPointsForSecondParam) * 1e-03;
tStart = tic;
LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
for i=1:NumberOfPointsForFirstParam
eval(sprintf('SidebandBeam.Detuning = %d;', DetuningArray(i)));
for j=1:NumberOfPointsForSecondParam
eval(sprintf('CoolingBeam.Waist = %d;', BeamWaistArray(j)));
eval(sprintf('SidebandBeam.Waist = %d;', BeamWaistArray(j)));
[LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOT2D.runSimulation(Oven);
end
end
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
if MOT2D.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
MOT2D.Results = LoadingRate;
SaveFolder = [MOT2D.SaveDirectory filesep 'Results'];
Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = MOT2D;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString = 'Sideband Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Beam Waist (mW)';
OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-10;
% OptionsStruct.ZLabelString = 'Enhancement Factor (\eta)';
OptionsStruct.ZLabelString = 'Loading rate (x 10^{10} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForTwoParameterScan(DetuningArray, BeamWaistArray, LoadingRateArray, options{:})
%%
MOT2D.NumberOfAtoms = 10000;
MOT2D.TotalPower = 0.4;
MOT2D.MagneticGradient = 0.4; 0;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power = 0.2;
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power = 0.2;
NumberOfPointsForFirstParam = 20; %iterations of the simulation
NumberOfPointsForSecondParam = 20;
DetuningArray = linspace(-0.5, -2.5, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth;
BeamWaistArray = linspace(10, 25, NumberOfPointsForSecondParam) * 1e-03;
tStart = tic;
LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
for i=1:NumberOfPointsForFirstParam
eval(sprintf('CoolingBeam.Detuning = %d;', DetuningArray(i)));
eval(sprintf('SidebandBeam.Detuning = %d;', DetuningArray(i) - (1.0 * Helper.PhysicsConstants.BlueLinewidth)));
for j=1:NumberOfPointsForSecondParam
eval(sprintf('CoolingBeam.Waist = %d;', BeamWaistArray(j)));
eval(sprintf('SidebandBeam.Waist = %d;', BeamWaistArray(j)));
[LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOT2D.runSimulation(Oven);
end
end
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
if MOT2D.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
MOT2D.Results = LoadingRate;
SaveFolder = [MOT2D.SaveDirectory filesep 'Results'];
Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = MOT2D;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString = 'Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Beam Waist (mW)';
OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-10;
% OptionsStruct.ZLabelString = 'Enhancement Factor (\eta)';
OptionsStruct.ZLabelString = 'Loading rate (x 10^{10} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForTwoParameterScan(DetuningArray, BeamWaistArray, LoadingRateArray, options{:})

33
MOT Simulator/+Scripts/scanForSidebandEnhancement.m Normal file
View File

@ -0,0 +1,33 @@
function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = scanForSidebandEnhancement(ovenObj, MOTobj, CBBeamName, SBBeamName, SBDetuningArray, SBPowerArray)
Beams = MOTobj.Beams;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, CBBeamName), Beams)};
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, SBBeamName), Beams)};
NumberOfPointsForFirstParam = length(SBDetuningArray);
NumberOfPointsForSecondParam = length(SBPowerArray);
LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
for i=1:NumberOfPointsForFirstParam
eval(sprintf('SidebandBeam.%s = %d;', 'Detuning', SBDetuningArray(i)));
for j=1:NumberOfPointsForSecondParam
eval(sprintf('SidebandBeam.%s = %d;', 'Power', SBPowerArray(j)));
eval(sprintf('CoolingBeam.%s = %d;', 'Power', MOTobj.TotalPower - SBPowerArray(j)));
[LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOTobj.runSimulation(ovenObj);
end
end
if MOTobj.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
LoadingRate.Errors = StandardErrorArray;
LoadingRate.CI = ConfidenceIntervalArray;
MOTobj.Results = LoadingRate;
SaveFolder = [MOTobj.SaveDirectory filesep 'Results'];
Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = MOTobj;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
end

27
MOT Simulator/+Simulator/+Scan/doOneParameter.m Normal file
View File

@ -0,0 +1,27 @@
function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOneParameter(ovenObj, MOTobj, BeamName, BeamParameter, ParameterArray)
Beams = MOTobj.Beams;
Beam = Beams{cellfun(@(x) strcmpi(x.Alias, BeamName), Beams)};
NumberOfPointsForParam = length(ParameterArray);
LoadingRateArray = zeros(1,NumberOfPointsForParam);
StandardErrorArray = zeros(1,NumberOfPointsForParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2);
for i=1:NumberOfPointsForParam
eval(sprintf('Beam.%s = %d;', BeamParameter, ParameterArray(i)));
[LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i,:)] = MOTobj.runSimulation(ovenObj);
end
if MOTobj.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
LoadingRate.Errors = StandardErrorArray;
LoadingRate.CI = ConfidenceIntervalArray;
MOTobj.Results = LoadingRate;
SaveFolder = [MOTobj.SaveDirectory filesep 'Results'];
Filename = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = MOTobj;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
end

22
MOT Simulator/+Simulator/+Scan/doThreeParameters.m Normal file
View File

@ -0,0 +1,22 @@
function LoadingRateArray = doThreeParameters(ovenObj, MOTobj, BeamName, FirstBeamParameter, FirstParameterArray, ...
SecondBeamParameter, SecondParameterArray, ThirdBeamParameter, ThirdParameterArray)
NumberOfPointsForThirdParam = length(ThirdParameterArray);
LoadingRateArray = {};
for i=1:NumberOfPointsForThirdParam
eval(sprintf('MOTobj.%s = %d;', ThirdBeamParameter, ThirdParameterArray(i)));
LoadingRateArray{end+1} = Simulator.Scan.doTwoParameters(ovenObj, MOTobj, BeamName, FirstBeamParameter, FirstParameterArray, SecondBeamParameter, SecondParameterArray);
end
if MOTobj.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
MOTobj.Results = LoadingRate;
SaveFolder = [MOTobj.SaveDirectory filesep 'Results'];
Filename = ['ThreeParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = MOTobj;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
end

31
MOT Simulator/+Simulator/+Scan/doTwoParameters.m Normal file
View File

@ -0,0 +1,31 @@
function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwoParameters(ovenObj, MOTobj, BeamName, FirstBeamParameter, ...
FirstParameterArray, SecondBeamParameter, SecondParameterArray)
Beams = MOTobj.Beams;
Beam = Beams{cellfun(@(x) strcmpi(x.Alias, BeamName), Beams)};
NumberOfPointsForFirstParam = length(FirstParameterArray);
NumberOfPointsForSecondParam = length(SecondParameterArray);
LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
StandardErrorArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
for i=1:NumberOfPointsForFirstParam
eval(sprintf('Beam.%s = %d;', FirstBeamParameter, FirstParameterArray(i)));
for j=1:NumberOfPointsForSecondParam
eval(sprintf('Beam.%s = %d;', SecondBeamParameter, SecondParameterArray(j)));
[LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOTobj.runSimulation(ovenObj);
end
end
if MOTobj.DoSave
LoadingRate = struct;
LoadingRate.Values = LoadingRateArray;
LoadingRate.Errors = StandardErrorArray;
LoadingRate.CI = ConfidenceIntervalArray;
MOTobj.Results = LoadingRate;
SaveFolder = [MOTobj.SaveDirectory filesep 'Results'];
Filename = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
eval([sprintf('%s_Object', Filename) ' = MOTobj;']);
mkdir(SaveFolder);
save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
end
end

203
MOT Simulator/+Simulator/@Beams/Beams.m Normal file
View File

@ -0,0 +1,203 @@
classdef Beams < handle & matlab.mixin.Copyable
properties (Access = private)
BlueBeamDefault = struct('Alias', 'Blue', ...
'Power', 400e-3, ...
'Detuning', -1.64*Helper.PhysicsConstants.BlueLinewidth, ...
'Radius', 17.5e-3, ...
'Waist', 15e-3, ...
'WaveNumber',2*pi/Helper.PhysicsConstants.BlueWavelength, ...
'Linewidth', Helper.PhysicsConstants.BlueLinewidth, ...
'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ...
((Helper.PhysicsConstants.PlanckConstantReduced * ...
Helper.PhysicsConstants.SpeedOfLight * ...
Helper.PhysicsConstants.BlueLinewidth) / (Helper.PhysicsConstants.BlueWavelength)^3));
BlueSidebandBeamDefault = struct('Alias', 'BlueSideband', ...
'Power', 400e-3, ...
'Detuning', -3*Helper.PhysicsConstants.BlueLinewidth, ...
'Radius', 17.5e-3, ...
'Waist', 15e-3, ...
'WaveNumber',2*pi/Helper.PhysicsConstants.BlueWavelength, ...
'Linewidth', Helper.PhysicsConstants.BlueLinewidth, ...
'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ...
((Helper.PhysicsConstants.PlanckConstantReduced * ...
Helper.PhysicsConstants.SpeedOfLight * ...
Helper.PhysicsConstants.BlueLinewidth) / (Helper.PhysicsConstants.BlueWavelength)^3));
PushBeamDefault = struct('Alias', 'Push', ...
'Power', 25e-3 , ...
'Detuning', 104.2*Helper.PhysicsConstants.RedLinewidth, ...
'Radius', 1.2e-03, ...
'Waist', 1.0e-03, ...
'WaveNumber',2*pi/Helper.PhysicsConstants.RedWavelength, ...
'Linewidth', Helper.PhysicsConstants.RedLinewidth, ...
'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ...
((Helper.PhysicsConstants.PlanckConstantReduced * ...
Helper.PhysicsConstants.SpeedOfLight * ...
Helper.PhysicsConstants.RedLinewidth) / (Helper.PhysicsConstants.RedWavelength)^3));
RedBeamDefault = struct('Alias', 'Red', ...
'Power', 100e-3 , ...
'Detuning', -1*Helper.PhysicsConstants.RedLinewidth, ...
'Radius', 1.2e-3, ...
'Waist', 12e-3 , ...
'WaveNumber',2*pi/Helper.PhysicsConstants.RedWavelength, ...
'Linewidth', Helper.PhysicsConstants.RedLinewidth, ...
'SaturationIntensity', 0.1 * (2 * pi^2 / 3) * ...
((Helper.PhysicsConstants.PlanckConstantReduced * ...
Helper.PhysicsConstants.SpeedOfLight * ...
Helper.PhysicsConstants.RedLinewidth) / (Helper.PhysicsConstants.RedWavelength)^3));
end
properties
Alias
Power;
Detuning;
Radius;
Waist;
WaveNumber;
SaturationIntensity;
Linewidth;
end
properties (Dependent)
SaturationParameter;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
methods
%% Class Constructor
function this = Beams(BeamName)
input = inputParser;
addRequired(input,'BeamName', @ischar);
parse(input, BeamName);
this.Alias = input.Results.BeamName;
switch this.Alias
case this.BlueBeamDefault.Alias
this.Power = this.BlueBeamDefault.Power;
this.Detuning = this.BlueBeamDefault.Detuning;
this.Radius = this.BlueBeamDefault.Radius;
this.Waist = this.BlueBeamDefault.Waist;
this.WaveNumber = this.BlueBeamDefault.WaveNumber;
this.Linewidth = this.BlueBeamDefault.Linewidth;
this.SaturationIntensity = this.BlueBeamDefault.SaturationIntensity;
case this.BlueSidebandBeamDefault.Alias
this.Power = this.BlueSidebandBeamDefault.Power;
this.Detuning = this.BlueSidebandBeamDefault.Detuning;
this.Radius = this.BlueSidebandBeamDefault.Radius;
this.Waist = this.BlueSidebandBeamDefault.Waist;
this.WaveNumber = this.BlueSidebandBeamDefault.WaveNumber;
this.Linewidth = this.BlueSidebandBeamDefault.Linewidth;
this.SaturationIntensity = this.BlueSidebandBeamDefault.SaturationIntensity;
case this.PushBeamDefault.Alias
this.Power = this.PushBeamDefault.Power;
this.Detuning = this.PushBeamDefault.Detuning;
this.Radius = this.PushBeamDefault.Radius;
this.Waist = this.PushBeamDefault.Waist;
this.WaveNumber = this.PushBeamDefault.WaveNumber;
this.Linewidth = this.PushBeamDefault.Linewidth;
this.SaturationIntensity = this.PushBeamDefault.SaturationIntensity;
case this.RedBeamDefault.Alias
this.Power = this.RedBeamDefault.Power;
this.Detuning = this.RedBeamDefault.Detuning;
this.Radius = this.RedBeamDefault.Radius;
this.Waist = this.RedBeamDefault.Waist;
this.WaveNumber = this.RedBeamDefault.WaveNumber;
this.Linewidth = this.RedBeamDefault.Linewidth;
this.SaturationIntensity = this.RedBeamDefault.SaturationIntensity;
otherwise
error('No such beam!')
end
end
end % - lifecycle
methods
function set.Power(this,val)
this.Power = val;
end
function ret = get.Power(this)
ret = this.Power;
end
function set.Detuning(this, val)
this.Detuning = val;
end
function ret = get.Detuning(this)
ret = this.Detuning;
end
function set.Radius(this, val)
this.Radius = val;
end
function ret = get.Radius(this)
ret = this.Radius;
end
function set.Waist(this, val)
this.Waist = val;
end
function ret = get.Waist(this)
ret = this.Waist;
end
function set.WaveNumber(this, val)
this.WaveNumber = val;
end
function ret = get.WaveNumber(this)
ret = this.WaveNumber;
end
function set.SaturationIntensity(this, val)
this.SaturationIntensity = val;
end
function ret = get.SaturationIntensity(this)
ret = this.SaturationIntensity;
end
function set.Linewidth(this, val)
this.Linewidth = val;
end
function ret = get.Linewidth(this)
ret = this.Linewidth;
end
end % - setters and getters
methods
function ret = get.SaturationParameter(this)
ret = 0.1 * (8 * this.Power) / (pi*this.Waist^2 * this.SaturationIntensity); % two beams are reflected
end
end % - getters for dependent properties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
methods(Access = protected)
function cp = copyElement(this)
% Shallow copy object
cp = copyElement@matlab.mixin.Copyable(this);
% Forces the setter to redefine the function handles to the new copied object
pl = properties(this);
for k = 1:length(pl)
sc = superclasses(this.(pl{k}));
if any(contains(sc,{'matlab.mixin.Copyable'}))
cp.(pl{k}) = this.(pl{k}).copy();
end
end
end
end
methods (Static)
% Creates an Instance of Class, ensures singleton behaviour (that there
% can only be one Instance of this class
function singleObj = getInstance(BeamName)
% Creates an Instance of Class, ensures singleton behaviour
persistent localObj;
if isempty(localObj) || ~isvalid(localObj)
localObj = Simulator.Beams(BeamName);
end
singleObj = localObj;
end
end
end

173
MOT Simulator/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m Normal file
View File

@ -0,0 +1,173 @@
classdef MOTCaptureProcess < handle & matlab.mixin.Copyable
properties (Access = public)
SimulationMode; % MOT type
TimeStep;
SimulationTime;
NumberOfAtoms;
ErrorEstimationMethod;
%Flags
SpontaneousEmission;
SidebandBeam;
PushBeam;
Gravity;
BackgroundCollision;
DebugMode;
DoSave;
SaveDirectory;
end
properties
Beams = {}; %Contains beam objects
end % - public
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
methods
%% Class Destructor (Clears object)
function this = MOTCaptureProcess(varargin)
p = inputParser;
p.KeepUnmatched = true;
addParameter(p, 'SimulationMode', '2D',...
@(x) any(strcmpi(x,{'2D','3D', 'Full'})));
addParameter(p, 'ErrorEstimationMethod', 'jackknife',...
@(x) any(strcmpi(x,{'jackknife','bootstrap'})));
addParameter(p, 'NumberOfAtoms', 5000,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'TimeStep', 10e-06,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'SimulationTime', 3e-03,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'SpontaneousEmission', false,...
@islogical);
addParameter(p, 'SidebandBeam', false,...
@islogical);
addParameter(p, 'PushBeam', false,...
@islogical);
addParameter(p, 'Gravity', false,...
@islogical);
addParameter(p, 'BackgroundCollision', false,...
@islogical);
addParameter(p, 'DebugMode', false,...
@islogical);
addParameter(p, 'SaveData', false,...
@islogical);
addParameter(p, 'SaveDirectory', pwd,...
@ischar);
p.parse(varargin{:});
this.SimulationMode = p.Results.SimulationMode;
this.ErrorEstimationMethod= p.Results.ErrorEstimationMethod;
this.NumberOfAtoms = p.Results.NumberOfAtoms;
this.TimeStep = p.Results.TimeStep;
this.SimulationTime = p.Results.SimulationTime;
this.SpontaneousEmission = p.Results.SpontaneousEmission;
this.SidebandBeam = p.Results.SidebandBeam;
this.PushBeam = p.Results.PushBeam;
this.Gravity = p.Results.Gravity;
this.BackgroundCollision = p.Results.BackgroundCollision;
this.DebugMode = p.Results.DebugMode;
this.DoSave = p.Results.SaveData;
this.SaveDirectory = p.Results.SaveDirectory;
switch this.SimulationMode
case "2D"
this.Beams{1} = Simulator.Beams('Blue');
this.Beams{2} = Simulator.Beams('BlueSideband');
this.Beams{3} = Simulator.Beams('Push');
case "3D"
this.Beams{1} = Simulator.Beams('Red');
% Development In progress
case "Full"
% Development In progress
end
end
end % - lifecycle
methods
function set.TimeStep(this, val)
assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!');
this.TimeStep = val;
end
function ret = get.TimeStep(this)
ret = this.TimeStep;
end
function set.SimulationTime(this, val)
% assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!');
this.SimulationTime = val;
end
function ret = get.SimulationTime(this)
ret = this.SimulationTime;
end
function set.NumberOfAtoms(this, val)
assert(val <= 50000, '!!Not time efficient to compute for atom numbers larger than 50,000!!');
this.NumberOfAtoms = val;
end
function ret = get.NumberOfAtoms(this)
ret = this.NumberOfAtoms;
end
function set.DebugMode(this, val)
this.DebugMode = val;
end
function ret = get.DebugMode(this)
ret = this.DebugMode;
end
function set.DoSave(this, val)
this.DoSave = val;
end
function ret = get.DoSave(this)
ret = this.DoSave;
end
function set.SaveDirectory(this, val)
this.SaveDirectory = val;
end
function ret = get.SaveDirectory(this)
ret = this.SaveDirectory;
end
end % - setters and getters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
methods(Access = protected)
function cp = copyElement(this)
% Shallow copy object
cp = copyElement@matlab.mixin.Copyable(this);
% Forces the setter to redefine the function handles to the new copied object
pl = properties(this);
for k = 1:length(pl)
sc = superclasses(this.(pl{k}));
if any(contains(sc,{'matlab.mixin.Copyable'}))
cp.(pl{k}) = this.(pl{k}).copy();
end
end
end
end
methods (Static)
% Creates an Instance of Class, ensures singleton behaviour (that there
% can only be one Instance of this class
function singleObj = getInstance(varargin)
% Creates an Instance of Class, ensures singleton behaviour
persistent localObj;
if isempty(localObj) || ~isvalid(localObj)
localObj = Simulator.MOTCaptureProcess(varargin{:});
end
singleObj = localObj;
end
end
end

177
MOT Simulator/+Simulator/@Oven/Oven.m Normal file
View File

@ -0,0 +1,177 @@
classdef Oven < Simulator.MOTCaptureProcess & matlab.mixin.Copyable
properties (Access = private)
OvenDefaults = struct('NozzleLength', 60e-3, ...
'NozzleRadius', 2.60e-3, ...
'OvenTemperature', 1000);
end
properties (Access = public)
NozzleLength;
NozzleRadius;
OvenTemperature;
VelocityCutoff;
ClausingFactor;
ReducedClausingFactor;
ReducedFlux;
end
properties (Dependent)
ExitDivergence;
Beta;
OvenDistance;
OvenTemperatureinKelvin;
AverageVelocity;
AtomicBeamDensity;
MeanFreePath;
CollisionTime;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
methods
function this = Oven(varargin)
this@Simulator.MOTCaptureProcess('SimulationMode', '2D', varargin{:});
this.NozzleLength = this.OvenDefaults.NozzleLength;
this.NozzleRadius = this.OvenDefaults.NozzleRadius;
this.OvenTemperature = this.OvenDefaults.OvenTemperature;
this.ClausingFactor = this.calculateClausingFactor();
[this.ReducedClausingFactor, ~] = this.calculateReducedClausingFactor();
end
function restoreDefaults(this)
this.NozzleLength = this.OvenDefaults.NozzleLength;
this.NozzleRadius = this.OvenDefaults.NozzleRadius;
this.OvenTemperature = this.OvenDefaults.OvenTemperature;
this.ClausingFactor = this.calculateClausingFactor();
[this.ReducedClausingFactor, ~] = this.calculateReducedClausingFactor();
end
end % - lifecycle
methods
function set.NozzleLength(this,val)
this.NozzleLength = val;
end
function ret = get.NozzleLength(this)
ret = this.NozzleLength;
end
function set.NozzleRadius(this,val)
this.NozzleRadius = val;
end
function ret = get.NozzleRadius(this)
ret = this.NozzleRadius;
end
function set.OvenTemperature(this,val)
this.OvenTemperature = val;
end
function ret = get.OvenTemperature(this)
ret = this.OvenTemperature;
end
function set.VelocityCutoff(this,val)
this.VelocityCutoff = val;
end
function ret = get.VelocityCutoff(this)
ret = this.VelocityCutoff;
end
function set.ClausingFactor(this,val)
this.ClausingFactor = val;
end
function ret = get.ClausingFactor(this)
ret = this.ClausingFactor;
end
function set.ReducedClausingFactor(this,val)
this.ReducedClausingFactor = val;
end
function ret = get.ReducedClausingFactor(this)
ret = this.ReducedClausingFactor;
end
function set.ReducedFlux(this,val)
this.ReducedFlux = val;
end
function ret = get.ReducedFlux(this)
ret = this.ReducedFlux;
end
end % - setters and getters
methods
function ret = get.Beta(this)
ret = 2 * (this.NozzleRadius/this.NozzleLength);
end
function ret = get.OvenDistance(this)
ApertureCut = max(2.5e-3,this.NozzleRadius);
ret = (25+12.5)*1e-3 + (this.NozzleRadius + ApertureCut) / tan(15/360 * 2 * pi);
end
function ret = get.ExitDivergence(this)
Theta_Nozzle = atan((this.NozzleRadius + 0.035/sqrt(2))/this.OvenDistance); % The angle of capture region towards the oven nozzle
Theta_Aperture = 15/360*2*pi; % The limitation angle of the second aperture in the oven
ret = min(Theta_Nozzle,Theta_Aperture);
end
function ret = get.OvenTemperatureinKelvin(this)
ret = this.OvenTemperature + Helper.PhysicsConstants.ZeroKelvin;
end
function ret = get.AverageVelocity(this)
%See Background collision probability section in Barbiero
ret = sqrt((8 * pi * Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass));
end
function ret = get.AtomicBeamDensity(this)
%See Background collision probability section in Barbiero
ret = this.calculateFreeMolecularRegimeFlux / (this.AverageVelocity * pi * (this.NozzleRadius)^2);
end
function ret = get.MeanFreePath(this)
% Cross section = pi ( 2 * Van-der-waals radius of Dy)^2;
% Van-der-waals radius of Dy = 281e-12
%See Expected atomic flux section and Background collision probability section in Barbiero
ret = 1/(sqrt(2) * (pi * (2*281e-12)^2) * this.AtomicBeamDensity);
end
function ret = get.CollisionTime(this)
ret = 3 * this.MeanFreePath/this.AverageVelocity; %See Background collision probability section in Barbiero
end
end % - getters for dependent properties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
methods(Access = protected)
function cp = copyElement(this)
% Shallow copy object
cp = copyElement@matlab.mixin.Copyable(this);
% Forces the setter to redefine the function handles to the new copied object
pl = properties(this);
for k = 1:length(pl)
sc = superclasses(this.(pl{k}));
if any(contains(sc,{'matlab.mixin.Copyable'}))
cp.(pl{k}) = this.(pl{k}).copy();
end
end
end
end
methods (Static)
% Creates an Instance of Class, ensures singleton behaviour (that there
% can only be one Instance of this class
function singleObj = getInstance(varargin)
% Creates an Instance of Class, ensures singleton behaviour
persistent localObj;
if isempty(localObj) || ~isvalid(localObj)
localObj = Simulator.Oven();
end
singleObj = localObj;
end
end
end

37
MOT Simulator/+Simulator/@Oven/angularDistributionFunction.m Normal file
View File

@ -0,0 +1,37 @@
function ret = angularDistributionFunction(this, theta)
%This function calculate the angle distribution of atoms coming out
%from a single channel.
KnudsenNumber = this.MeanFreePath/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))));
eta0 = alpha;
eta1 = 1 - alpha;
delta = (eta0./sqrt(2*KnudsenNumber*(eta1-eta0)))./sqrt(cos(theta));
F = 2/sqrt(pi)* (1-eta1)/eta0 * delta.* exp( -(delta*eta1/eta0).^2 );
q = this.Beta^-1 * tan(theta);
R = acos(q) - (q .* sqrt(1-q.^2));
if abs(q) >= 1
t = linspace(0,1,10000);
S = sum(sqrt(1-t.^2).* ( erf(delta.*(1 + (t.*(eta1-eta0)./(q.*eta0)) ))-erf(delta)))*(t(2)-t(1));
if S == 0 || isnan(S)
ret = eta0*cos(theta);
else
ret = eta0*cos(theta)+ 2/sqrt(pi)*eta0*cos(theta) * (exp(delta.^2)/delta) * S;
end
else
t = linspace(0,q,10000);
S = sum(sqrt(1-t.^2).* ( erf(delta.*(1 + (t.*(eta1-eta0)./(q.*eta0)) ))-erf(delta)))*(t(2)-t(1));
if isnan(S)
S=0;
end
ret = 2/sqrt(pi)*eta0*cos(theta)*(exp(delta.^2)/delta) * (R./2*(erf(delta*eta1/eta0)-erf(delta)+F)+S)+eta0*cos(theta);
end
end

9
MOT Simulator/+Simulator/@Oven/calculateClausingFactor.m Normal file
View File

@ -0,0 +1,9 @@
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 = ClausingFactorApproximation;
end

11
MOT Simulator/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m Normal file
View File

@ -0,0 +1,11 @@
function ret = calculateFreeMolecularRegimeFlux(this)
%This function calculate the total flux of atoms coming out from a tube
%See Expected atomic flux section in Barbiero
Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature
Dy164DensityinOven = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin);
ret = 1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2;
% Removed the Helper.PhysicsConstants.Dy164IsotopicAbundance multiplication
% 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.
end

17
MOT Simulator/+Simulator/@Oven/calculateReducedClausingFactor.m Normal file
View File

@ -0,0 +1,17 @@
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.ExitDivergence
ReducedClausingFactor = ReducedClausingFactor + (2 * pi * sin(ThetaArray(p)) * AngularDistribution(p) * (ThetaArray(2)-ThetaArray(1)));
end
end
ReducedClausingFactor = ReducedClausingFactor / pi;
end

7
MOT Simulator/+Simulator/@Oven/initialPositionSampling.m Normal file
View File

@ -0,0 +1,7 @@
function ret = initialPositionSampling(this)
n = this.NumberOfAtoms;
phi = 2 * pi * rand(n,1);
rho = this.Beta * 0.5 * this.NozzleLength * sqrt(rand(n,1));
ret = [-this.OvenDistance * ones(n,1), rho.*cos(phi), rho.*sin(phi)];
end

60
MOT Simulator/+Simulator/@Oven/initialVelocitySampling.m Normal file
View File

@ -0,0 +1,60 @@
function ret = initialVelocitySampling(this, MOTObj)
n = this.NumberOfAtoms;
% Calculate Calculate Capture velocity --> Introduce velocity cutoff
MOTObj.CaptureVelocity = MOTObj.calculateCaptureVelocity(this, [-this.OvenDistance,0,0], [1,0,0]);
this.VelocityCutoff = 1.05 * MOTObj.CaptureVelocity(1); % Should be the magnitude of the 3-D velocity vector but since here the obtained capture
% velocity is only along the x-axis, we take the first term which is the x-component of the velocity.
[ReducedClausingFactor, NormalizationConstantForAngularDistribution] = this.calculateReducedClausingFactor();
this.ReducedClausingFactor = ReducedClausingFactor;
VelocityDistribution = @(velocity) sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin))^3 ...
* velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* this.OvenTemperatureinKelvin)));
c = integral(VelocityDistribution, 0, this.VelocityCutoff) / integral(VelocityDistribution, 0, inf);
this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
ret = zeros(n,3);
SampledVelocityMagnitude = zeros(n,1);
SampledPolarAngle = zeros(n,1);
SampledAzimuthalAngle = zeros(n,1);
MostProbableVelocity = sqrt((3 * Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperature) / Helper.PhysicsConstants.Dy164Mass); % For v * f(v) distribution
if MostProbableVelocity > this.VelocityCutoff
MaximumVelocityAllowed = this.VelocityCutoff;
else
MaximumVelocityAllowed = MostProbableVelocity;
end
NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed);
parfor i = 1:n
% Rejection Sampling of speed
y = rand(1);
x = this.VelocityCutoff * rand(1);
while y > ((NormalizationConstantForVelocityDistribution)^-1 * this.velocityDistributionFunction(x)) %As long as this loop condition is satisfied, reject the corresponding x value
y = rand(1);
x = this.VelocityCutoff * rand(1);
end
SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample
% Rejection Sampling of polar angle
w = rand(1);
z = this.ExitDivergence * rand(1);
while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value
w = rand(1);
z = this.ExitDivergence * rand(1);
end
SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
% Sampling of azimuthal angle
SampledAzimuthalAngle(i)= 2 * pi * rand(1);
ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ...
SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))];
end
end

5
MOT Simulator/+Simulator/@Oven/velocityDistributionFunction.m Normal file
View File

@ -0,0 +1,5 @@
function ret = velocityDistributionFunction(this, velocity)
ret = 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)));
end

191
MOT Simulator/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m Normal file
View File

@ -0,0 +1,191 @@
classdef TwoDimensionalMOT < Simulator.MOTCaptureProcess & matlab.mixin.Copyable
properties (Access = private)
MagneticGradienDefault = 0.40; % T/m
ExitDivergenceDefault = 15e-3;
DistanceBetweenPushBeamAnd3DMOTCenterDefault = 0;
PushBeamDistanceDefault = 0.32;
end
properties (Access = public)
TotalPower;
LandegFactor;
MagneticSubLevel;
MagneticGradient;
CaptureVelocity;
ExitDivergence;
DistanceBetweenPushBeamAnd3DMOTCenter;
PushBeamDistance;
TimeSpentInInteractionRegion;
ParticleDynamicalQuantities;
InitialParameters;
BootstrapSampleLength;
BootstrapSampleNumber;
Results;
end
methods
function this = TwoDimensionalMOT(varargin)
this@Simulator.MOTCaptureProcess('SimulationMode', '2D', varargin{:});
p = inputParser;
p.KeepUnmatched = true;
addParameter(p, 'TotalPower', 0.8,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'LandegFactor', 1,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'MagneticSubLevel', 1,...
@(x) assert(isnumeric(x) && isscalar(x)));
addParameter(p, 'MagneticGradient', 0.38,...
@(x) assert(isnumeric(x) && isscalar(x)));
p.parse(varargin{:});
this.TotalPower = p.Results.TotalPower;
this.LandegFactor = p.Results.LandegFactor;
this.MagneticSubLevel = p.Results.MagneticSubLevel;
this.MagneticGradient = p.Results.MagneticGradient;
this.ExitDivergence = this.ExitDivergenceDefault;
this.DistanceBetweenPushBeamAnd3DMOTCenter = this.DistanceBetweenPushBeamAnd3DMOTCenterDefault;
this.PushBeamDistance = this.PushBeamDistanceDefault;
this.InitialParameters = struct;
this.InitialParameters.TotalPower = this.TotalPower;
this.InitialParameters.LandegFactor = this.LandegFactor;
this.InitialParameters.MagneticSubLevel= this.MagneticSubLevel;
this.InitialParameters.MagneticGradient= this.MagneticGradient;
this.BootstrapSampleLength = 0.5 * this.NumberOfAtoms;
this.BootstrapSampleNumber = 1000;
end
function restoreDefaults(this)
this.TotalPower = this.InitialParameters.TotalPower;
this.LandegFactor = this.InitialParameters.LandegFactor;
this.MagneticSubLevel = this.InitialParameters.MagneticSubLevel;
this.MagneticGradient = this.InitialParameters.MagneticGradient;
this.ExitDivergence = this.ExitDivergenceDefault;
this.DistanceBetweenPushBeamAnd3DMOTCenter = this.DistanceBetweenPushBeamAnd3DMOTCenterDefault;
this.PushBeamDistance = this.PushBeamDistanceDefault;
end
end % - lifecycle
methods
function set.TotalPower(this,val)
this.TotalPower = val;
end
function ret = get.TotalPower(this)
ret = this.TotalPower;
end
function set.LandegFactor(this,val)
this.LandegFactor = val;
end
function ret = get.LandegFactor(this)
ret = this.LandegFactor;
end
function set.MagneticSubLevel(this,val)
this.MagneticSubLevel = val;
end
function ret = get.MagneticSubLevel(this)
ret = this.MagneticSubLevel;
end
function set.MagneticGradient(this,val)
this.MagneticGradient = val;
end
function ret = get.MagneticGradient(this)
ret = this.MagneticGradient;
end
function set.CaptureVelocity(this,val)
this.CaptureVelocity = val;
end
function ret = get.CaptureVelocity(this)
ret = this.CaptureVelocity;
end
function set.ExitDivergence(this,val)
this.ExitDivergence = val;
end
function ret = get.ExitDivergence(this)
ret = this.ExitDivergence;
end
function set.DistanceBetweenPushBeamAnd3DMOTCenter(this,val)
this.DistanceBetweenPushBeamAnd3DMOTCenter = val;
end
function ret = get.DistanceBetweenPushBeamAnd3DMOTCenter(this)
ret = this.DistanceBetweenPushBeamAnd3DMOTCenter;
end
function set.PushBeamDistance(this,val)
this.PushBeamDistance = val;
end
function ret = get.PushBeamDistance(this)
ret = this.PushBeamDistance;
end
function set.TimeSpentInInteractionRegion(this,val)
this.TimeSpentInInteractionRegion = val;
end
function ret = get.TimeSpentInInteractionRegion(this)
ret = this.TimeSpentInInteractionRegion;
end
function set.ParticleDynamicalQuantities(this,val)
this.ParticleDynamicalQuantities = val;
end
function ret = get.ParticleDynamicalQuantities(this)
ret = this.ParticleDynamicalQuantities;
end
function set.InitialParameters(this,val)
this.InitialParameters = val;
end
function ret = get.InitialParameters(this)
ret = this.InitialParameters;
end
function set.BootstrapSampleLength(this,val)
this.BootstrapSampleLength = val;
end
function ret = get.BootstrapSampleLength(this)
ret = this.BootstrapSampleLength;
end
function set.BootstrapSampleNumber(this,val)
this.BootstrapSampleNumber = val;
end
function ret = get.BootstrapSampleNumber(this)
ret = this.BootstrapSampleNumber;
end
function set.Results(this, val)
this.Results = val;
end
function ret = get.Results(this)
ret = this.Results;
end
end % - setters and getters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
methods(Access = protected)
function cp = copyElement(this)
% Shallow copy object
cp = copyElement@matlab.mixin.Copyable(this);
% Forces the setter to redefine the function handles to the new copied object
pl = properties(this);
for k = 1:length(pl)
sc = superclasses(this.(pl{k}));
if any(contains(sc,{'matlab.mixin.Copyable'}))
cp.(pl{k}) = this.(pl{k}).copy();
end
end
end
end
methods (Static)
% Creates an Instance of Class, ensures singleton behaviour (that there
% can only be one Instance of this class
function singleObj = getInstance(varargin)
% Creates an Instance of Class, ensures singleton behaviour
persistent localObj;
if isempty(localObj) || ~isvalid(localObj)
localObj = Simulator.TwoDimensionalMOT(varargin{:});
end
singleObj = localObj;
end
end
end

35
MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m Normal file
View File

@ -0,0 +1,35 @@
function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
% is the distance between the chamber center and the cross point of push beam and z-axis (along the gravity)
WaveVectorEndPoint = [0, 1, this.DistanceBetweenPushBeamAnd3DMOTCenter/this.PushBeamDistance];
WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint);
Origin=[0,0,0];
PushBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), this.Beams)};
PushBeamDetuning = PushBeamObj.Detuning;
PushBeamRadius = PushBeamObj.Radius;
PushBeamWaist = PushBeamObj.Waist;
PushBeamWaveNumber = PushBeamObj.WaveNumber;
PushBeamLinewidth = PushBeamObj.Linewidth;
PushBeamSaturationParameter = 0.25 * PushBeamObj.SaturationParameter;
SaturationIntensity = this.calculateLocalSaturationIntensity(PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, PushBeamRadius, PushBeamWaist);
DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * PushBeamWaveNumber;
Detuning = PushBeamDetuning - DopplerShift;
s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./PushBeamLinewidth).^2));
a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(PushBeamLinewidth * 0.5);
a_push = a_sat .* (s_push/(1+s_push));
if this.SpontaneousEmission
a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, PushBeamLinewidth, PushBeamWaveNumber);
else
a_scatter = [0,0,0];
end
a_total = a_push + a_scatter;
ret = a_total(1:3);
end

15
MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m Normal file
View File

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

22
MOT Simulator/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m Normal file
View File

@ -0,0 +1,22 @@
function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstimation(this, ovenObj, NumberOfLoadedAtoms)
n = this.NumberOfAtoms;
SampleLength = this.BootstrapSampleLength;
NumberOfBootsrapSamples = this.BootstrapSampleNumber;
MeanCaptureRatioInEachSample = zeros(1,NumberOfBootsrapSamples);
for SampleNumber = 1:NumberOfBootsrapSamples
BoostrapSample = datasample(NumberOfLoadedAtoms, SampleLength); % Sample with replacement
MeanCaptureRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
end
LoadingRate = mean(MeanCaptureRatioInEachSample) * ovenObj.ReducedFlux;
Variance = 0; % Bootstrap Estimate of Variance
for SampleNumber = 1:NumberOfBootsrapSamples
Variance = Variance + (MeanCaptureRatioInEachSample(SampleNumber) - mean(MeanCaptureRatioInEachSample))^2;
end
StandardError = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * ovenObj.ReducedFlux;
ts = tinv([0.025 0.975],NumberOfBootsrapSamples-1); % T-Score
ConfidenceInterval = LoadingRate + ts*StandardError; % 95% Confidence Intervals
end

22
MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m Normal file
View File

@ -0,0 +1,22 @@
function ret = calculateCaptureVelocity(this, ovenObj, PositionVector, VelocityVector)
VelocityUnitVector = VelocityVector./norm(VelocityVector);
UpperLimit = 500;
LowerLimit = 0;
for Index = 1:500
InitialVelocity = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector;
ParticleDynamicalQuantities = this.solver(PositionVector, InitialVelocity);
FinalPositionVector = ParticleDynamicalQuantities(end, 1:3);
if rssq(FinalPositionVector) <= ovenObj.OvenDistance
LowerLimit = 0.5 * (UpperLimit + LowerLimit);
else
UpperLimit = 0.5 * (UpperLimit + LowerLimit);
end
if UpperLimit - LowerLimit < 1
ret = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector;
break;
end
end
end

53
MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m Normal file
View File

@ -0,0 +1,53 @@
function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ovenObj)
n = this.NumberOfAtoms;
DynamicalQuantities = this.ParticleDynamicalQuantities;
CollisionEvents = zeros(1, n);
NumberOfLoadedAtoms = zeros(1, n);
% Include the stochastic process of background collisions
for AtomIndex = 1:n
this.TimeSpentInInteractionRegion(AtomIndex) = this.computeTimeSpentInInteractionRegion(squeeze(DynamicalQuantities(AtomIndex,:,1:3)));
CollisionEvents(AtomIndex) = this.computeCollisionProbability(ovenObj, this.TimeSpentInInteractionRegion(AtomIndex));
end
% Count the number of loaded atoms subject to conditions
switch this.ErrorEstimationMethod
case 'bootstrap'
NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep);
NumberOfLoadedAtoms = zeros(1, NumberOfTimeSteps);
LoadedAtomIndices = [];
for TimeIndex = 1:NumberOfTimeSteps
if TimeIndex ~= 1
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
end
for AtomIndex = 1:n
Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
if this.exitCondition(Position, CollisionEvents(AtomIndex))
if ~ismember(AtomIndex, LoadedAtomIndices)
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
LoadedAtomIndices(end+1) = AtomIndex;
end
else
if ismember(AtomIndex, LoadedAtomIndices)
NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) - 1;
LoadedAtomIndices(LoadedAtomIndices==AtomIndex) = [];
end
end
end
end
[LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(ovenObj, NumberOfLoadedAtoms);
case 'jackknife'
for AtomIndex = 1:n
if AtomIndex ~= 1
NumberOfLoadedAtoms(AtomIndex) = NumberOfLoadedAtoms(AtomIndex-1);
end
Position = squeeze(DynamicalQuantities(AtomIndex, end, 1:3))';
if this.exitCondition(Position, CollisionEvents(AtomIndex))
NumberOfLoadedAtoms(AtomIndex) = NumberOfLoadedAtoms(AtomIndex) + 1;
end
end
[LoadingRate, StandardError, ConfidenceInterval] = jackknifeErrorEstimation(this, ovenObj, NumberOfLoadedAtoms);
end
end

10
MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m Normal file
View File

@ -0,0 +1,10 @@
function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVector, WaveVectorOrigin, WaveVectorEndPoint, BeamRadius, BeamWaist)
DistanceBetweenAtomAndLaserBeamAxis = Helper.calculateDistanceFromPointToLine(PositionVector, WaveVectorOrigin, WaveVectorEndPoint);
if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius
ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2);
else
ret = 0;
end
end

104
MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m Normal file
View File

@ -0,0 +1,104 @@
function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
CoolingBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), this.Beams)};
CoolingBeamDetuning = CoolingBeamObj.Detuning;
CoolingBeamRadius = CoolingBeamObj.Radius;
CoolingBeamWaist = CoolingBeamObj.Waist;
CoolingBeamWaveNumber = CoolingBeamObj.WaveNumber;
CoolingBeamLinewidth = CoolingBeamObj.Linewidth;
CoolingBeamSaturationParameter = CoolingBeamObj.SaturationParameter;
SidebandBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), this.Beams)};
SidebandDetuning = SidebandBeamObj.Detuning;
SidebandBeamRadius = SidebandBeamObj.Radius;
SidebandBeamWaist = SidebandBeamObj.Waist;
SidebandSaturationParameter = SidebandBeamObj.SaturationParameter;
WaveVectorEndPoint = zeros(2,3);
WaveVectorEndPoint(1,:) = [1,0,1];
WaveVectorEndPoint(1,:) = WaveVectorEndPoint(1,1:3)/norm(WaveVectorEndPoint(1,:));
WaveVectorEndPoint(2,:) = [-1,0,1];
WaveVectorEndPoint(2,:) = WaveVectorEndPoint(2,1:3)/norm(WaveVectorEndPoint(2,:));
Sigma = [1,-1];
Origin = [0,0,0];
% Calculate the Saturation Intensity at the specified point along its Gaussian Profile
CoolingBeamLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(0.25 * CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), CoolingBeamRadius, CoolingBeamWaist), ...
this.calculateLocalSaturationIntensity(0.25 * CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), CoolingBeamRadius, CoolingBeamWaist)];
SidebandLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(0.25 * SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), SidebandBeamRadius, SidebandBeamWaist), ...
this.calculateLocalSaturationIntensity(0.25 * SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), SidebandBeamRadius, SidebandBeamWaist)];
TotalAcceleration = zeros(1,3);
Delta_Cooling = [0,0,0,0];
Delta_Sideband = [0,0,0,0];
for i = 1:2
LocalMagneticField = this.magneticFieldForMOT(PositionVector);
B = sign(dot(LocalMagneticField(1:3), WaveVectorEndPoint(i,:))) * LocalMagneticField(4);
ZeemanShift = this.LandegFactor * this.MagneticSubLevel * (Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced) * B;
DopplerShift = dot(WaveVectorEndPoint(i,:), VelocityVector) * CoolingBeamWaveNumber;
Delta_Cooling(i*2-1) = CoolingBeamDetuning + DopplerShift + (ZeemanShift * Sigma(i));
Delta_Cooling(i*2) = CoolingBeamDetuning - DopplerShift - (ZeemanShift * Sigma(i));
if this.SidebandBeam
Delta_Sideband(i*2-1) = SidebandDetuning + DopplerShift + (ZeemanShift * Sigma(i));
Delta_Sideband(i*2) = SidebandDetuning - DopplerShift - (ZeemanShift * Sigma(i));
end
end
SaturationParameter = [0,0,0,0,0,0,0,0];
for i = 1:2
SaturationParameter(2*i-1) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/CoolingBeamLinewidth)^2);
SaturationParameter(2*i) = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i) /CoolingBeamLinewidth)^2);
if this.SidebandBeam
SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/CoolingBeamLinewidth)^2);
SaturationParameter(2*i+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i)/CoolingBeamLinewidth)^2);
end
end
TotalSaturationParameter = sum(SaturationParameter);
for i = 1:2
a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(CoolingBeamLinewidth * 0.5);
a_1 = a_sat .* (SaturationParameter(2*i-1)/(1 + TotalSaturationParameter));
a_2 = a_sat .* (SaturationParameter(2*i) /(1 + TotalSaturationParameter));
if this.SpontaneousEmission
a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber) + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber);
else
a_scattering = [0,0,0];
end
if this.SidebandBeam
a_1 = a_1 + a_sat .* (SaturationParameter(2*i-1+4)/(1 + TotalSaturationParameter));
a_2 = a_2 + a_sat .* (SaturationParameter(2*i+4) /(1 + TotalSaturationParameter));
if this.SpontaneousEmission
a_scattering = a_scattering + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber) + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber);
else
a_scattering = [0,0,0];
end
end
TotalAcceleration = TotalAcceleration + (a_2 - a_1) + a_scattering;
end
if this.PushBeam
TotalAcceleration = TotalAcceleration + this.accelerationDueToPushBeam(PositionVector, VelocityVector);
end
ret = TotalAcceleration(1:3);
end

9
MOT Simulator/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m Normal file
View File

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

20
MOT Simulator/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m Normal file
View File

@ -0,0 +1,20 @@
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;
CoolingBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), this.Beams)};
NumberOfTimeSteps = int64(this.SimulationTime/this.TimeStep);
for n = 1:(NumberOfTimeSteps - 1)
dr = Helper.calculateDistanceFromPointToLine(r(n+1, :), [0 0 0], [0 0 1]);
if dr < CoolingBeamObj.Radius
A = 1;
else
A = 0;
end
T = T + A * this.TimeStep;
end
end

10
MOT Simulator/+Simulator/@TwoDimensionalMOT/exitCondition.m Normal file
View File

@ -0,0 +1,10 @@
function ret = exitCondition(this, PositionVector, CollisionEvent)
d = Helper.calculateDistanceFromPointToLine(PositionVector, [0 0 0], [0 1 0]);
y = PositionVector(2);
DivergenceAngle = atan(d/abs(y));
if (y >= 0) && (DivergenceAngle <= this.ExitDivergence) && ~CollisionEvent
ret = true;
else
ret = false;
end
end

50
MOT Simulator/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m Normal file
View File

@ -0,0 +1,50 @@
function [LoadingRate, StandardError, ConfidenceInterval] = jackknifeErrorEstimation(this, ovenObj, NumberOfLoadedAtoms)
n = this.NumberOfAtoms;
Autocorrelation = zeros(1, n);
for i = 1:n-1
FirstTerm = 0;
SecondTerm = 0;
for j = 1:n-i
FirstTerm = FirstTerm + NumberOfLoadedAtoms(j) / j;
SecondTerm = SecondTerm + (NumberOfLoadedAtoms(i+j)) / (i+j);
Autocorrelation(i) = Autocorrelation(i) + ((NumberOfLoadedAtoms(j) / j) .*(NumberOfLoadedAtoms(i+j) / (i+j)));
end
Autocorrelation(i) = (1/(n-i)) * (Autocorrelation(i) - ((1/(n-i)) * FirstTerm * SecondTerm));
end
if Autocorrelation(1)~=0
Autocorrelation = Autocorrelation./Autocorrelation(1);
x = linspace(1,n,n);
[FitParams,~] = fit(x',Autocorrelation',"exp(-x/tau)", 'Startpoint', 100);
CorrelationFactor = FitParams.tau;
SampleLength = 2*CorrelationFactor+1;
NumberOfJackknifeSamples = floor(n/SampleLength);
CaptureRatioInEachSample = zeros(1,NumberOfJackknifeSamples);
SampleNumberLimit = min(NumberOfJackknifeSamples-1,5);
for i=1:NumberOfJackknifeSamples-SampleNumberLimit
CaptureRatioInEachSample(i) = NumberOfLoadedAtoms(n-ceil((i-1)*SampleLength))/(n-ceil((i-1)*SampleLength));
end
MeanCaptureRatio = sum(CaptureRatioInEachSample) / (NumberOfJackknifeSamples-SampleNumberLimit);
LoadingRate = MeanCaptureRatio * ovenObj.ReducedFlux;
Variance=0;
for i=1:NumberOfJackknifeSamples-SampleNumberLimit
Variance=Variance+(CaptureRatioInEachSample(i) - MeanCaptureRatio)^2;
end
StandardError = sqrt(Variance/(NumberOfJackknifeSamples-SampleNumberLimit));
ConfidenceInterval = LoadingRate + 1.96*StandardError; % 95% Confidence Intervals
else
LoadingRate = nan;
StandardError = nan;
ConfidenceInterval = [nan nan];
end
end

8
MOT Simulator/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m Normal file
View File

@ -0,0 +1,8 @@
function ret = magneticFieldForMOT(this, r)
ret = zeros(1,4);
alpha = this.MagneticGradient;
ret(1) = r(3)*alpha;
ret(2) = 0;
ret(3) = r(1)*alpha;
ret(4) = sqrt(ret(1)^2+ret(2)^2+ret(3)^2);
end

25
MOT Simulator/+Simulator/@TwoDimensionalMOT/runSimulation.m Normal file
View File

@ -0,0 +1,25 @@
function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this, ovenObj)
if this.NumberOfAtoms ~= ovenObj.NumberOfAtoms
ovenObj.NumberOfAtoms = this.NumberOfAtoms;
end
%% - Sampling for initial positions and velocities
% - sampling the position distribution
Positions = ovenObj.initialPositionSampling();
% - sampling the velocity distribution
Velocities = ovenObj.initialVelocitySampling(this);
%% Solve ODE
progressbar = Helper.parforNotifications();
progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating 2-D MOT capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']);
% calculate the final position of the atoms
DynamicalQuantities = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
parfor Index = 1:this.NumberOfAtoms
DynamicalQuantities(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(ovenObj);
end

39
MOT Simulator/+Simulator/@TwoDimensionalMOT/solver.m Normal file
View File

@ -0,0 +1,39 @@
function ParticleDynamicalQuantities = 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);
for i=1:int64(this.SimulationTime/this.TimeStep)
ParticleDynamicalQuantities(i,1:3) = InitialPosition;
ParticleDynamicalQuantities(i,4:6) = InitialVelocity;
rt = InitialPosition;
vt = InitialVelocity;
ga1 = this.calculateTotalAcceleration(rt,vt) + g;
gv1 = vt .* this.TimeStep;
rt = rt + 0.5 * gv1;
vt = vt + 0.5 * ga1 .* this.TimeStep;
ga2 = this.calculateTotalAcceleration(rt,vt) + g;
gv2 = vt .* this.TimeStep;
rt = rt + 0.5 * gv2;
vt = vt + 0.5 * ga2 .* this.TimeStep;
ga3 = this.calculateTotalAcceleration(rt,vt) + g;
gv3 = vt .* this.TimeStep;
rt = rt + 0.5 * gv3;
vt = vt + ga3 .* this.TimeStep;
ga4 = this.calculateTotalAcceleration(rt,vt) + g;
gv4 = vt .* this.TimeStep;
InitialPosition = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6;
InitialVelocity = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6;
end
end

427
MOT Simulator/test_MOTCaptureProcessSimulation.m Normal file
View File

@ -0,0 +1,427 @@
%% This script is testing the functionalities of the MOT Capture Process Simulation Classes
%
% Important: Run only sectionwise!!
%% - Testing the MOTCaptureProcess-Class
% - Create MOTCaptureProcess object with specified options
% - Automatically creates Beams objects
OptionsStruct = struct;
OptionsStruct.ErrorEstimationMethod = 'bootstrap'; % 'jackknife' | 'bootstrap'
OptionsStruct.NumberOfAtoms = 5000;
OptionsStruct.TimeStep = 50e-06; % in s
OptionsStruct.SimulationTime = 4e-03; % in s
OptionsStruct.SpontaneousEmission = true;
OptionsStruct.SidebandBeam = true;
OptionsStruct.PushBeam = true;
OptionsStruct.Gravity = true;
OptionsStruct.BackgroundCollision = true;
OptionsStruct.SaveData = true;
% OptionsStruct.SaveDirectory = '';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
Oven = Simulator.Oven(options{:});
MOT2D = Simulator.TwoDimensionalMOT(options{:});
Beams = MOT2D.Beams;
%% - Run Simulation
MOT2D.NumberOfAtoms = 10000;
MOT2D.SidebandBeam = true;
MOT2D.PushBeam = false;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power = 0.2;
CoolingBeam.Waist = 20e-03;
CoolingBeam.Detuning = -1.33*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power = 0.2;
SidebandBeam.Waist = 20e-03;
SidebandBeam.Detuning = -2.66*Helper.PhysicsConstants.BlueLinewidth;
PushBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), Beams)};
PushBeam.Power = 0.025;
PushBeam.Waist = 0.81e-03;
PushBeam.Detuning = 0;
[LoadingRate, ~] = MOT2D.runSimulation(Oven);
%% - Plot initial distribution
% - sampling the position distribution
InitialPositions = Oven.initialPositionSampling();
% - sampling the velocity distribution
InitialVelocities = Oven.initialVelocitySampling(MOT2D);
NumberOfBins = 100;
Plotter.plotPositionAndVelocitySampling(NumberOfBins, InitialPositions, InitialVelocities);
%% - Plot distributions of magnitude and direction of initial velocities
NumberOfBins = 50;
Plotter.plotInitialVeloctiySamplingVsAngle(Oven, MOT2D, NumberOfBins)
%% - Plot Magnetic Field
XAxisRange = [-5 5];
YAxisRange = [-5 5];
ZAxisRange = [-5 5];
Plotter.visualizeMagneticField(MOT2D, XAxisRange, YAxisRange, ZAxisRange)
%% - Plot MFP & VP for different temperatures
TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius
Plotter.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
%% - Plot the Free Molecular Flux for different temperatures
Temperature = [950, 1000, 1050]; % Temperature
Plotter.plotFreeMolecularFluxVsTemp(Oven,Temperature)
%% - Plot Angular Distribution for different Beta
Beta = [0.5, 0.1 , 0.05, 0.02, 0.01]; %Beta = 2 * radius / length of the tube
Plotter.plotAngularDistributionForDifferentBeta(Oven, Beta)
%% - Plot Capture Velocity
Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot!
%% - Plot Phase Space with Acceleration Field
MOT2D.SidebandBeam = false;
MOT2D.MagneticGradient = 0.4;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power = 0.3;
CoolingBeam.Detuning = -1.64*Helper.PhysicsConstants.BlueLinewidth;
CoolingBeam.Waist = 15e-03;
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power = 0.5;
SidebandBeam.Detuning = -4*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam.Waist = 15e-03;
MOT2D.NumberOfAtoms = 50;
MinimumVelocity = 0;
MaximumVelocity = 150;
NumberOfBins = 200; %Along each axis
IncidentAtomDirection = 0*2*pi/360;
IncidentAtomPosition = 0;
Plotter.plotPhaseSpaceWithAccelerationField(Oven, MOT2D, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
%% - Plot Trajectories along the 3 directions
MOT2D.NumberOfAtoms = 100;
MOT2D.MagneticGradient = 0.42;
MaximumVelocity = 150;
IncidentAtomDirection = 0*2*pi/360;
IncidentAtomPosition = 0;
%% - Positions
Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, 'PlotPositions', true);
%% - Velocities
Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, 'PlotVelocities', true);
%% - Scan parameters: One-Parameter Scan
MOT2D.NumberOfAtoms = 5000;
MOT2D.TotalPower = 0.4;
MOT2D.SidebandBeam = false;
MOT2D.PushBeam = false;
NumberOfPointsForFirstParam = 5; %iterations of the simulation
ParameterArray = linspace(0.1, 1.0, NumberOfPointsForFirstParam) * MOT2D.TotalPower;
tStart = tic;
[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.Scan.doOneParameter(Oven, MOT2D, 'Blue', 'Power', ParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForParameter = 1000;
OptionsStruct.XLabelString = 'Cooling Beam Power (mW)';
OptionsStruct.RescalingFactorForYQuantity = 1e-10;
OptionsStruct.ErrorsForYQuantity = true;
OptionsStruct.ErrorsArray = StandardErrorArray;
OptionsStruct.CIForYQuantity = true;
OptionsStruct.CIArray = ConfidenceIntervalArray;
OptionsStruct.RemoveOutliers = true;
OptionsStruct.YLabelString = 'Loading rate (x 10^{10} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForOneParameterScan(ParameterArray, LoadingRateArray, options{:})
clear OptionsStruct
%% - Scan parameters: One-Parameter Scan
MOT2D.NumberOfAtoms = 10000;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power = 0.4;
MOT2D.SidebandBeam = false;
MOT2D.PushBeam = false;
% ParameterArray = [10 20 30 40 50 60 70 80 90 100];
ParameterArray = [500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500];
NumberOfPointsForParam = length(ParameterArray); %iterations of the simulation
LoadingRateArray = zeros(1,NumberOfPointsForParam);
StandardErrorArray = zeros(1,NumberOfPointsForParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2);
tStart = tic;
for i=1:NumberOfPointsForParam
MOT2D.BootstrapSampleLength = ParameterArray(i);
[LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i,:)] = MOT2D.runSimulation(Oven);
end
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForParameter = 1;
OptionsStruct.XLabelString = 'Bootstrap Sample Length';
OptionsStruct.RescalingFactorForYQuantity = 1e-10;
OptionsStruct.ErrorsForYQuantity = true;
OptionsStruct.ErrorsArray = StandardErrorArray;
OptionsStruct.CIForYQuantity = true;
OptionsStruct.CIArray = ConfidenceIntervalArray;
OptionsStruct.RemoveOutliers = false;
OptionsStruct.YLabelString = 'Loading rate (x 10^{10} atoms/s)';
OptionsStruct.TitleString = sprintf('Cooling Beam Power = %d (mW); Magnetic Gradient = %.0f (G/cm)', CoolingBeam.Power*1000, MOT2D.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForOneParameterScan(ParameterArray, LoadingRateArray, options{:})
MeanLR = mean(LoadingRateArray(:)) * 1e-10;
yline(MeanLR, 'LineStyle', '--', 'Linewidth', 2.5)
textstring = [sprintf('%1.2e', MeanLR * 1e+10) ' atoms'];
% txt = text((ParameterArray(2) + 0.05*ParameterArray(2)), (max(MeanLR) + 0.05*MeanLR), textstring, 'Interpreter','latex', 'FontSize', 14);
% xlim([0 100])
ylim([0 3.5])
clear OptionsStruct
%% - Scan parameters: Two-Parameter Scan
% COOLING BEAM POWER VS DETUNING
MOT2D.NumberOfAtoms = 5000;
MOT2D.TotalPower = 0.6;
NumberOfPointsForFirstParam = 10; %iterations of the simulation
NumberOfPointsForSecondParam = 10;
FirstParameterArray = linspace(-0.5, -2.5, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth;
SecondParameterArray = linspace(0.3, 1.0, NumberOfPointsForSecondParam) * MOT2D.TotalPower;
tStart = tic;
[LoadingRateArray, ~, ~] = Simulator.Scan.doTwoParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, 'Power', SecondParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString = 'Cooling Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Cooling Beam Power (mW)';
OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-9;
OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)';
OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, LoadingRateArray, options{:})
clear OptionsStruct
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% COOLING BEAM WAIST VS DETUNING
MOT2D.NumberOfAtoms = 20000;
MOT2D.MagneticGradient = 0.40;
MOT2D.SidebandBeam = false;
MOT2D.PushBeam = false;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power = 0.4;
NumberOfPointsForFirstParam = 10; %iterations of the simulation
NumberOfPointsForSecondParam = 10;
FirstParameterArray = linspace(-0.5, -2.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth;
SecondParameterArray = linspace(10, 25, NumberOfPointsForSecondParam) * 1e-03;
tStart = tic;
[LoadingRateArray, ~, ~] = Simulator.Scan.doTwoParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, 'Waist', SecondParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString = 'Cooling Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Cooling Beam Waist (mm)';
OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-9;
OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)';
OptionsStruct.TitleString = sprintf('Cooling Beam Power = %d (mW); Magnetic Gradient = %.0f (G/cm)', CoolingBeam.Power*1000, MOT2D.MagneticGradient * 100);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, LoadingRateArray, options{:})
clear OptionsStruct
%% - Scan parameters: Three-Parameter Scan
% COOLING BEAM WAIST VS DETUNING FOR DIFFERENT MAGNETIC FIELD GRADIENTS
MOT2D.NumberOfAtoms = 10000;
MOT2D.SidebandBeam = false;
MOT2D.PushBeam = false;
MOT2D.BackgroundCollision = false;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power = 0.4;
NumberOfPointsForFirstParam = 10; %iterations of the simulation
NumberOfPointsForSecondParam = 10;
NumberOfPointsForThirdParam = 6;
FirstParameterArray = linspace(-0.5, -2.0, NumberOfPointsForFirstParam) * Helper.PhysicsConstants.BlueLinewidth;
SecondParameterArray = linspace(10, 25, NumberOfPointsForSecondParam) * 1e-03;
ThirdParameterArray = linspace(30, 50, NumberOfPointsForThirdParam) * 1e-02;
MOT2D.BootstrapSampleLength = 500;
tStart = tic;
LoadingRateArray = Simulator.Scan.doThreeParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, ...
'Waist', SecondParameterArray, ...
'MagneticGradient', ThirdParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
% - Plot results
OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString = 'Cooling Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter = 1000;
OptionsStruct.YLabelString = 'Cooling Beam Waist (mm)';
OptionsStruct.RescalingFactorForThirdParameter = 100;
OptionsStruct.RescalingFactorForQuantityOfInterest = 1e-9;
OptionsStruct.ZLabelString = 'Loading rate (x 10^{9} atoms/s)';
OptionsStruct.PlotTitleString = 'Magnetic Gradient = %.0f (G/cm)';
OptionsStruct.FigureTitleString = sprintf('Oven-2DMOT Distance = %.1f (mm); Cooling Beam Power = %d (mW)', Oven.OvenDistance * 1000, CoolingBeam.Power*1000);
options = Helper.convertstruct2cell(OptionsStruct);
Plotter.plotResultForThreeParameterScan(FirstParameterArray, SecondParameterArray, ThirdParameterArray, LoadingRateArray, options{:})
clear OptionsStruct
%% Local Saturation Intensity distribution
WaveVectorEndPoint = zeros(2,3);
WaveVectorEndPoint(1,:) = [1,0,1];
WaveVectorEndPoint(1,:) = WaveVectorEndPoint(1,1:3)/norm(WaveVectorEndPoint(1,:));
WaveVectorEndPoint(2,:) = [-1,0,1];
WaveVectorEndPoint(2,:) = WaveVectorEndPoint(2,1:3)/norm(WaveVectorEndPoint(2,:));
Origin = [0,0,0];
BeamNumber = 2; %Selects one of the two wave vectors defined above
BeamRadius = 17.5e-03;
BeamWaist = 15e-03;
Power = 0.4;
CoolingBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
SaturationIntensity = CoolingBeam.SaturationIntensity;
SaturationParameter = 0.1 * (8 * Power) / (pi*BeamWaist^2 * SaturationIntensity); % two beams are reflected
n = 10000;
xmin = -BeamRadius;
xmax = BeamRadius;
x = xmin+rand(n,1)*(xmax-xmin);
y = ones(n,1) * 0;
zmin = -BeamRadius;
zmax = BeamRadius;
z = zmin+rand(n,1)*(zmax-zmin);
% t = 2*pi*rand(n,1);
% r = BeamRadius*sqrt(rand(n,1));
% x = r.*cos(t);
% y = ones(n,1) * 0;
% z = r.*sin(t);
PositionVector = horzcat(x, y, z); %scatter3(zeros(n,1), y, z)
CoolingBeamLocalSaturationIntensity = @(x) MOT2D.calculateLocalSaturationIntensity(0.25 * SaturationParameter, x, Origin, WaveVectorEndPoint(BeamNumber,:), BeamRadius, BeamWaist);
IntensityProfile = zeros(n,1);
for i=1:n
IntensityProfile(i) = CoolingBeamLocalSaturationIntensity(PositionVector(i, :));
end
v = IntensityProfile; % Extract intensity value
rows = 35;
columns = 35;
Image = zeros(rows, columns);
for k = 1 : length(x)
row = ceil((x(k) - min(x)) * columns / (max(x) - min(x)));
column = ceil((z(k) - min(z)) * rows / (max(z) - min(z)));
if (row == 0)
row = 1;
end
if (column == 0)
column = 1;
end
Image(row, column) = v(k);
end
f_h = Helper.getFigureByTag('Intensity Profile');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Intensity Profile';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
imagesc(linspace(min(x),max(x),row) * 1e+03, linspace(min(z),max(z),column) * 1e+03, Image);
set(gca,'YDir','normal');
hXLabel = xlabel('x-direction (distance in mm)');
hYLabel = ylabel('z-direction (distance in mm)');
shading flat;
c = colorbar;
c.Label.String= 'Local Saturation Intensity';
c.Label.FontSize = 14;
hTitle = sgtitle('Intensity Distribution');
set([hXLabel, hYLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();
%% Beam Shape in Three dimensions
f_h = Helper.getFigureByTag('Intensity Profile');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children'))
clf(f_h);
end
f_h.Name = 'Intensity Profile';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
[xq,zq] = meshgrid(linspace(-BeamRadius, BeamRadius, n), linspace(-BeamRadius, BeamRadius, n));
vq = griddata(x,z,v,xq,zq);
mesh(xq,zq,vq)
hold on
plot3(x,z,v,'o', 'MarkerSize', 1.5)
hXLabel = xlabel('x-direction (distance in mm)');
hYLabel = ylabel('z-direction (distance in mm)');
shading flat;
c = colorbar;
c.Label.String= 'Local Saturation Intensity';
c.Label.FontSize = 14;
hTitle = sgtitle('Intensity Distribution');
set([hXLabel, hYLabel] , ...
'FontSize' , 14 );
set( hTitle , ...
'FontSize' , 18 );
Helper.bringFiguresWithTagInForeground();