From 33390f8230d6ba30ccfeced58bb2f148be496156 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 12 Jun 2024 19:09:15 +0200 Subject: [PATCH] Newly renamed, restructured Simulator code. --- .../+Helper/PhysicsConstants.m | 46 + .../+Helper/parforNotifications.m | 148 ++++ Dipolar Gas Simulator/+Plotter/makeMovie.m | 77 ++ Dipolar Gas Simulator/+Plotter/plotLive.m | 47 + .../+Plotter/visualizeSpace.m | 94 ++ Dipolar Gas Simulator/+Scripts/analyze.m | 50 ++ .../@Calculator/calculateChemicalPotential.m | 28 + .../@Calculator/calculateEnergyComponents.m | 35 + .../calculateNormalizedResiduals.m | 24 + .../calculateNumericalHankelTransform.m | 39 + .../@Calculator/calculateOrderParameter.m | 58 ++ .../@Calculator/calculatePhaseCoherence.m | 18 + .../@Calculator/calculateTotalEnergy.m | 31 + .../@Calculator/calculateVDCutoff.m | 33 + .../+Simulator/@DipolarGas/initialize.m | 22 + .../+Simulator/@DipolarGas/runSimulation.m | 35 + .../+Simulator/@DipolarGas/setupParameters.m | 101 +++ .../+Simulator/@DipolarGas/setupSpace.m | 33 + .../+Simulator/@DipolarGas/setupSpaceRadial.m | 307 +++++++ .../@DipolarGas/setupWavefunction.m | 23 + .../splitStepFourier.m | 90 ++ MOT Simulator/+Helper/ImageSelection.class | Bin 0 -> 1092 bytes MOT Simulator/+Helper/ImageSelection.java | 38 + MOT Simulator/+Helper/PhysicsConstants.m | 46 + .../+Helper/bringFiguresWithTagInForeground.m | 15 + .../calculateDistanceFromPointToLine.m | 10 + MOT Simulator/+Helper/convertstruct2cell.m | 6 + MOT Simulator/+Helper/findAllZeroCrossings.m | 18 + MOT Simulator/+Helper/getFigureByTag.m | 191 ++++ MOT Simulator/+Helper/ode5.m | 92 ++ MOT Simulator/+Helper/onenoteccdata.m | 55 ++ MOT Simulator/+Helper/parforNotifications.m | 148 ++++ MOT Simulator/+Helper/screencapture.m | 820 ++++++++++++++++++ .../plotAngularDistributionForDifferentBeta.m | 73 ++ .../+Plotter/plotCaptureVelocityVsAngle.m | 37 + .../+Plotter/plotConfidenceIntervalRegion.m | 18 + .../+Plotter/plotDynamicalQuantities.m | 85 ++ .../+Plotter/plotFreeMolecularFluxVsTemp.m | 54 ++ .../plotInitialVeloctiySamplingVsAngle.m | 72 ++ .../plotMeanFreePathAndVapourPressureVsTemp.m | 43 + .../plotPhaseSpaceWithAccelerationField.m | 84 ++ .../plotPositionAndVelocitySampling.m | 55 ++ .../+Plotter/plotResultForOneParameterScan.m | 90 ++ .../plotResultForThreeParameterScan.m | 82 ++ .../+Plotter/plotResultForTwoParameterScan.m | 74 ++ .../+Plotter/visualizeMagneticField.m | 72 ++ .../optimizingForSidebandEnhancement.m | 236 +++++ .../+Scripts/scanForSidebandEnhancement.m | 33 + .../+Simulator/+Scan/doOneParameter.m | 27 + .../+Simulator/+Scan/doThreeParameters.m | 22 + .../+Simulator/+Scan/doTwoParameters.m | 31 + MOT Simulator/+Simulator/@Beams/Beams.m | 203 +++++ .../@MOTCaptureProcess/MOTCaptureProcess.m | 173 ++++ MOT Simulator/+Simulator/@Oven/Oven.m | 177 ++++ .../@Oven/angularDistributionFunction.m | 37 + .../@Oven/calculateClausingFactor.m | 9 + .../@Oven/calculateFreeMolecularRegimeFlux.m | 11 + .../@Oven/calculateReducedClausingFactor.m | 17 + .../@Oven/initialPositionSampling.m | 7 + .../@Oven/initialVelocitySampling.m | 60 ++ .../@Oven/velocityDistributionFunction.m | 5 + .../@TwoDimensionalMOT/TwoDimensionalMOT.m | 191 ++++ .../accelerationDueToPushBeam.m | 35 + ...elerationDueToSpontaneousEmissionProcess.m | 15 + .../bootstrapErrorEstimation.m | 22 + .../calculateCaptureVelocity.m | 22 + .../@TwoDimensionalMOT/calculateLoadingRate.m | 53 ++ .../calculateLocalSaturationIntensity.m | 10 + .../calculateTotalAcceleration.m | 104 +++ .../computeCollisionProbability.m | 9 + .../computeTimeSpentInInteractionRegion.m | 20 + .../@TwoDimensionalMOT/exitCondition.m | 10 + .../jackknifeErrorEstimation.m | 50 ++ .../@TwoDimensionalMOT/magneticFieldForMOT.m | 8 + .../@TwoDimensionalMOT/runSimulation.m | 25 + .../+Simulator/@TwoDimensionalMOT/solver.m | 39 + .../test_MOTCaptureProcessSimulation.m | 427 +++++++++ 77 files changed, 5705 insertions(+) create mode 100644 Dipolar Gas Simulator/+Helper/PhysicsConstants.m create mode 100644 Dipolar Gas Simulator/+Helper/parforNotifications.m create mode 100644 Dipolar Gas Simulator/+Plotter/makeMovie.m create mode 100644 Dipolar Gas Simulator/+Plotter/plotLive.m create mode 100644 Dipolar Gas Simulator/+Plotter/visualizeSpace.m create mode 100644 Dipolar Gas Simulator/+Scripts/analyze.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m create mode 100644 Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m create mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m create mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m create mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m create mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m create mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m create mode 100644 Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m create mode 100644 Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m create mode 100644 MOT Simulator/+Helper/ImageSelection.class create mode 100644 MOT Simulator/+Helper/ImageSelection.java create mode 100644 MOT Simulator/+Helper/PhysicsConstants.m create mode 100644 MOT Simulator/+Helper/bringFiguresWithTagInForeground.m create mode 100644 MOT Simulator/+Helper/calculateDistanceFromPointToLine.m create mode 100644 MOT Simulator/+Helper/convertstruct2cell.m create mode 100644 MOT Simulator/+Helper/findAllZeroCrossings.m create mode 100644 MOT Simulator/+Helper/getFigureByTag.m create mode 100644 MOT Simulator/+Helper/ode5.m create mode 100644 MOT Simulator/+Helper/onenoteccdata.m create mode 100644 MOT Simulator/+Helper/parforNotifications.m create mode 100644 MOT Simulator/+Helper/screencapture.m create mode 100644 MOT Simulator/+Plotter/plotAngularDistributionForDifferentBeta.m create mode 100644 MOT Simulator/+Plotter/plotCaptureVelocityVsAngle.m create mode 100644 MOT Simulator/+Plotter/plotConfidenceIntervalRegion.m create mode 100644 MOT Simulator/+Plotter/plotDynamicalQuantities.m create mode 100644 MOT Simulator/+Plotter/plotFreeMolecularFluxVsTemp.m create mode 100644 MOT Simulator/+Plotter/plotInitialVeloctiySamplingVsAngle.m create mode 100644 MOT Simulator/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m create mode 100644 MOT Simulator/+Plotter/plotPhaseSpaceWithAccelerationField.m create mode 100644 MOT Simulator/+Plotter/plotPositionAndVelocitySampling.m create mode 100644 MOT Simulator/+Plotter/plotResultForOneParameterScan.m create mode 100644 MOT Simulator/+Plotter/plotResultForThreeParameterScan.m create mode 100644 MOT Simulator/+Plotter/plotResultForTwoParameterScan.m create mode 100644 MOT Simulator/+Plotter/visualizeMagneticField.m create mode 100644 MOT Simulator/+Scripts/optimizingForSidebandEnhancement.m create mode 100644 MOT Simulator/+Scripts/scanForSidebandEnhancement.m create mode 100644 MOT Simulator/+Simulator/+Scan/doOneParameter.m create mode 100644 MOT Simulator/+Simulator/+Scan/doThreeParameters.m create mode 100644 MOT Simulator/+Simulator/+Scan/doTwoParameters.m create mode 100644 MOT Simulator/+Simulator/@Beams/Beams.m create mode 100644 MOT Simulator/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m create mode 100644 MOT Simulator/+Simulator/@Oven/Oven.m create mode 100644 MOT Simulator/+Simulator/@Oven/angularDistributionFunction.m create mode 100644 MOT Simulator/+Simulator/@Oven/calculateClausingFactor.m create mode 100644 MOT Simulator/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m create mode 100644 MOT Simulator/+Simulator/@Oven/calculateReducedClausingFactor.m create mode 100644 MOT Simulator/+Simulator/@Oven/initialPositionSampling.m create mode 100644 MOT Simulator/+Simulator/@Oven/initialVelocitySampling.m create mode 100644 MOT Simulator/+Simulator/@Oven/velocityDistributionFunction.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/exitCondition.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/runSimulation.m create mode 100644 MOT Simulator/+Simulator/@TwoDimensionalMOT/solver.m create mode 100644 MOT Simulator/test_MOTCaptureProcessSimulation.m diff --git a/Dipolar Gas Simulator/+Helper/PhysicsConstants.m b/Dipolar Gas Simulator/+Helper/PhysicsConstants.m new file mode 100644 index 0000000..f62dd0a --- /dev/null +++ b/Dipolar Gas Simulator/+Helper/PhysicsConstants.m @@ -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 diff --git a/Dipolar Gas Simulator/+Helper/parforNotifications.m b/Dipolar Gas Simulator/+Helper/parforNotifications.m new file mode 100644 index 0000000..4ad3af4 --- /dev/null +++ b/Dipolar Gas Simulator/+Helper/parforNotifications.m @@ -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 + + diff --git a/Dipolar Gas Simulator/+Plotter/makeMovie.m b/Dipolar Gas Simulator/+Plotter/makeMovie.m new file mode 100644 index 0000000..8842cbd --- /dev/null +++ b/Dipolar Gas Simulator/+Plotter/makeMovie.m @@ -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 \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Plotter/plotLive.m b/Dipolar Gas Simulator/+Plotter/plotLive.m new file mode 100644 index 0000000..1e3abb3 --- /dev/null +++ b/Dipolar Gas Simulator/+Plotter/plotLive.m @@ -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]); diff --git a/Dipolar Gas Simulator/+Plotter/visualizeSpace.m b/Dipolar Gas Simulator/+Plotter/visualizeSpace.m new file mode 100644 index 0000000..ac5646c --- /dev/null +++ b/Dipolar Gas Simulator/+Plotter/visualizeSpace.m @@ -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 \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Scripts/analyze.m b/Dipolar Gas Simulator/+Scripts/analyze.m new file mode 100644 index 0000000..df85ab4 --- /dev/null +++ b/Dipolar Gas Simulator/+Scripts/analyze.m @@ -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); diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m new file mode 100644 index 0000000..697ab38 --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m @@ -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; \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m new file mode 100644 index 0000000..49159ed --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m @@ -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) \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m new file mode 100644 index 0000000..66e5316 --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m @@ -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); \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m new file mode 100644 index 0000000..5e5be07 --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m @@ -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 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 \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m new file mode 100644 index 0000000..d970966 --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m @@ -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 \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m new file mode 100644 index 0000000..e0a5d3c --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m @@ -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 \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m new file mode 100644 index 0000000..79ade2a --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m @@ -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 \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m b/Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m new file mode 100644 index 0000000..c69ceb5 --- /dev/null +++ b/Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Helper/ImageSelection.class b/MOT Simulator/+Helper/ImageSelection.class new file mode 100644 index 0000000000000000000000000000000000000000..cfac41d4d552a6c83b184c77190049bc9ccc5c9a GIT binary patch literal 1092 zcma)5%Wl&^6g^`nabnywl(w{?v`zXT4`_e|3#eG|sDxC(A_b|ia*|2p%5{R{6uyOT zU;*V(2?XrjthEC>baV7;V=+==yzul^5H4_JLiqUj?<69oT_yd;PZbYY%wY3cKzIB%OV` zBLx=Y<}g#cH)yk2wjQZE8&jK(=LB~J3Z;LymY)eE?sr=xo!oXj`FOD3kp7O{a8;%w zgPmg`N{7I$5xUc4mZOQT?R9ET8hf%CP>}iXbyM~Nr|WUq*}r(B{a9ElmCxkEjMI;O zsSkR+t{=#j!pGa5D(|^Kdb8;s8>E+%1!lcF@SAeWQEOiaU93x&(kXaDJ&c7M7A#~j zX~DvTg$nWl*T=uvQ?GxbDOzo~yrQWJERV;P}Zs_nC9HJrT!tW|l&l&#*p}-H`1d-5?S03>np((?7CYjISJmVB^73MXbX5|Q? zQvC$&J#X}#F$3JKmam}YhwGwfEl+pH;EzIq5PO`xw0EqI z^2|}kJn@$>%SwVZrQ{;!7!~6JPoXL#jIpUOx5PNlO`^`?@oaNA`|WU6)W3=}=O{+S fyv}LrmrZ 0 + output = figure_handles; +end + +end \ No newline at end of file diff --git a/MOT Simulator/+Helper/calculateDistanceFromPointToLine.m b/MOT Simulator/+Helper/calculateDistanceFromPointToLine.m new file mode 100644 index 0000000..df5c8c6 --- /dev/null +++ b/MOT Simulator/+Helper/calculateDistanceFromPointToLine.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Helper/convertstruct2cell.m b/MOT Simulator/+Helper/convertstruct2cell.m new file mode 100644 index 0000000..90fdf2c --- /dev/null +++ b/MOT Simulator/+Helper/convertstruct2cell.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Helper/findAllZeroCrossings.m b/MOT Simulator/+Helper/findAllZeroCrossings.m new file mode 100644 index 0000000..4b8d9db --- /dev/null +++ b/MOT Simulator/+Helper/findAllZeroCrossings.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Helper/getFigureByTag.m b/MOT Simulator/+Helper/getFigureByTag.m new file mode 100644 index 0000000..8fc6bf9 --- /dev/null +++ b/MOT Simulator/+Helper/getFigureByTag.m @@ -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 + + + diff --git a/MOT Simulator/+Helper/ode5.m b/MOT Simulator/+Helper/ode5.m new file mode 100644 index 0000000..3eb003f --- /dev/null +++ b/MOT Simulator/+Helper/ode5.m @@ -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.'; diff --git a/MOT Simulator/+Helper/onenoteccdata.m b/MOT Simulator/+Helper/onenoteccdata.m new file mode 100644 index 0000000..5f3c177 --- /dev/null +++ b/MOT Simulator/+Helper/onenoteccdata.m @@ -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'); \ No newline at end of file diff --git a/MOT Simulator/+Helper/parforNotifications.m b/MOT Simulator/+Helper/parforNotifications.m new file mode 100644 index 0000000..4ad3af4 --- /dev/null +++ b/MOT Simulator/+Helper/parforNotifications.m @@ -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 + + diff --git a/MOT Simulator/+Helper/screencapture.m b/MOT Simulator/+Helper/screencapture.m new file mode 100644 index 0000000..206312e --- /dev/null +++ b/MOT Simulator/+Helper/screencapture.m @@ -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 MathWorks File Exchange + +% 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 + 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 + 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 + 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 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 + 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 \ No newline at end of file diff --git a/MOT Simulator/+Plotter/plotAngularDistributionForDifferentBeta.m b/MOT Simulator/+Plotter/plotAngularDistributionForDifferentBeta.m new file mode 100644 index 0000000..3812724 --- /dev/null +++ b/MOT Simulator/+Plotter/plotAngularDistributionForDifferentBeta.m @@ -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 diff --git a/MOT Simulator/+Plotter/plotCaptureVelocityVsAngle.m b/MOT Simulator/+Plotter/plotCaptureVelocityVsAngle.m new file mode 100644 index 0000000..92c9306 --- /dev/null +++ b/MOT Simulator/+Plotter/plotCaptureVelocityVsAngle.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Plotter/plotConfidenceIntervalRegion.m b/MOT Simulator/+Plotter/plotConfidenceIntervalRegion.m new file mode 100644 index 0000000..2769798 --- /dev/null +++ b/MOT Simulator/+Plotter/plotConfidenceIntervalRegion.m @@ -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 + \ No newline at end of file diff --git a/MOT Simulator/+Plotter/plotDynamicalQuantities.m b/MOT Simulator/+Plotter/plotDynamicalQuantities.m new file mode 100644 index 0000000..dad74a4 --- /dev/null +++ b/MOT Simulator/+Plotter/plotDynamicalQuantities.m @@ -0,0 +1,85 @@ +function plotDynamicalQuantities(OvenObj, MOTObj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, varargin) + + p = inputParser; + p.KeepUnmatched = true; + + addRequired(p, 'OvenObject' , @isobject); + addRequired(p, 'MOTObject' , @isobject); + addRequired(p, 'MaximumVelocity' , @(x) assert(isnumeric(x) && isscalar(x))); + addRequired(p, 'IncidentAtomDirection' , @(x) assert(isnumeric(x) && isscalar(x))); + addRequired(p, 'IncidentAtomPosition' , @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'PlotPositions' , false, @islogical); + addParameter(p, 'PlotVelocities' , false, @islogical); + + p.parse(OvenObj, MOTObj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, varargin{:}) + + MaximumVelocity = p.Results.MaximumVelocity; + IncidentAtomDirection = p.Results.IncidentAtomDirection; + IncidentAtomPosition = p.Results.IncidentAtomDirection; + PlotPositions = p.Results.PlotPositions; + PlotVelocities = p.Results.PlotVelocities; + + f_h = Helper.getFigureByTag('Trajectories Plot'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Trajectories Plot'; + f_h.Units = 'pixels'; + + set(0,'units','pixels'); + screensize = get(0,'ScreenSize'); + f_h.Position = [[screensize(3)/50 screensize(4)/3.5] 1870 600]; + + N = MOTObj.NumberOfAtoms; + Theta = IncidentAtomDirection; + z = IncidentAtomPosition; + + L = OvenObj.OvenDistance * 2; + T = MOTObj.SimulationTime; + tau = MOTObj.TimeStep; + + Y = linspace(0,MaximumVelocity,N); + DynamicalQuantities = zeros(length(Y),int64(T/tau),6); + for i=1:length(Y) + x =-L/2; + vx = Y(i)*cos(Theta); + vz = Y(i)*sin(Theta); + r = [x,0,z]; + v = [vx,0,vz]; + DynamicalQuantities(i,:,:) = MOTObj.solver(r, v); + end + + LabelStringArrayForPositions = {'x (mm)', 'y (mm)', 'z (mm)'}; + LabelStringArrayForVelocities = {'v_x (m/s)', 'v_y (m/s)', 'v_z (m/s)'}; + colors = { [0, 0.4470, 0.7410], [0.8500, 0.3250, 0.0980], [0.4940, 0.1840, 0.5560]}; + t = linspace(0, T, T/tau) * 1e+3; + for i=1:length(Y) + for j = 1:3 + if PlotPositions + subplot(1, 3, j) + hold on + plot(t, DynamicalQuantities(i,:,j) * 1e+3,'Color', colors{j},'linewidth',1.3) + hXLabel = xlabel('Time (ms)'); + hYLabel = ylabel(LabelStringArrayForPositions{j}); + hold off + set([hXLabel, hYLabel], 'FontSize', 14); + elseif PlotVelocities + subplot(1, 3, j) + hold on + plot(t, DynamicalQuantities(i,:,j+3),'Color', colors{j},'linewidth',1.3) + hXLabel = xlabel('Time (ms)'); + hYLabel = ylabel(LabelStringArrayForVelocities{j}); + hold off + set([hXLabel, hYLabel], 'FontSize', 14); + end + end + end + + hold off + hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", MOTObj.MagneticGradient)); + set(hTitle, 'FontSize', 18); + Helper.bringFiguresWithTagInForeground(); +end + diff --git a/MOT Simulator/+Plotter/plotFreeMolecularFluxVsTemp.m b/MOT Simulator/+Plotter/plotFreeMolecularFluxVsTemp.m new file mode 100644 index 0000000..46e5734 --- /dev/null +++ b/MOT Simulator/+Plotter/plotFreeMolecularFluxVsTemp.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Plotter/plotInitialVeloctiySamplingVsAngle.m b/MOT Simulator/+Plotter/plotInitialVeloctiySamplingVsAngle.m new file mode 100644 index 0000000..47b9479 --- /dev/null +++ b/MOT Simulator/+Plotter/plotInitialVeloctiySamplingVsAngle.m @@ -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 diff --git a/MOT Simulator/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m b/MOT Simulator/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m new file mode 100644 index 0000000..e7b9551 --- /dev/null +++ b/MOT Simulator/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m @@ -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 diff --git a/MOT Simulator/+Plotter/plotPhaseSpaceWithAccelerationField.m b/MOT Simulator/+Plotter/plotPhaseSpaceWithAccelerationField.m new file mode 100644 index 0000000..f198b01 --- /dev/null +++ b/MOT Simulator/+Plotter/plotPhaseSpaceWithAccelerationField.m @@ -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 + diff --git a/MOT Simulator/+Plotter/plotPositionAndVelocitySampling.m b/MOT Simulator/+Plotter/plotPositionAndVelocitySampling.m new file mode 100644 index 0000000..ad87281 --- /dev/null +++ b/MOT Simulator/+Plotter/plotPositionAndVelocitySampling.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Plotter/plotResultForOneParameterScan.m b/MOT Simulator/+Plotter/plotResultForOneParameterScan.m new file mode 100644 index 0000000..6ff5dcd --- /dev/null +++ b/MOT Simulator/+Plotter/plotResultForOneParameterScan.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Plotter/plotResultForThreeParameterScan.m b/MOT Simulator/+Plotter/plotResultForThreeParameterScan.m new file mode 100644 index 0000000..6867c48 --- /dev/null +++ b/MOT Simulator/+Plotter/plotResultForThreeParameterScan.m @@ -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', '', @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 \ No newline at end of file diff --git a/MOT Simulator/+Plotter/plotResultForTwoParameterScan.m b/MOT Simulator/+Plotter/plotResultForTwoParameterScan.m new file mode 100644 index 0000000..36141c0 --- /dev/null +++ b/MOT Simulator/+Plotter/plotResultForTwoParameterScan.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Plotter/visualizeMagneticField.m b/MOT Simulator/+Plotter/visualizeMagneticField.m new file mode 100644 index 0000000..2af9d23 --- /dev/null +++ b/MOT Simulator/+Plotter/visualizeMagneticField.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Scripts/optimizingForSidebandEnhancement.m b/MOT Simulator/+Scripts/optimizingForSidebandEnhancement.m new file mode 100644 index 0000000..4af4332 --- /dev/null +++ b/MOT Simulator/+Scripts/optimizingForSidebandEnhancement.m @@ -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{:}) \ No newline at end of file diff --git a/MOT Simulator/+Scripts/scanForSidebandEnhancement.m b/MOT Simulator/+Scripts/scanForSidebandEnhancement.m new file mode 100644 index 0000000..2b13fb7 --- /dev/null +++ b/MOT Simulator/+Scripts/scanForSidebandEnhancement.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/+Scan/doOneParameter.m b/MOT Simulator/+Simulator/+Scan/doOneParameter.m new file mode 100644 index 0000000..5ac8b56 --- /dev/null +++ b/MOT Simulator/+Simulator/+Scan/doOneParameter.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/+Scan/doThreeParameters.m b/MOT Simulator/+Simulator/+Scan/doThreeParameters.m new file mode 100644 index 0000000..4afdcf3 --- /dev/null +++ b/MOT Simulator/+Simulator/+Scan/doThreeParameters.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/+Scan/doTwoParameters.m b/MOT Simulator/+Simulator/+Scan/doTwoParameters.m new file mode 100644 index 0000000..092ca13 --- /dev/null +++ b/MOT Simulator/+Simulator/+Scan/doTwoParameters.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@Beams/Beams.m b/MOT Simulator/+Simulator/@Beams/Beams.m new file mode 100644 index 0000000..7889f1d --- /dev/null +++ b/MOT Simulator/+Simulator/@Beams/Beams.m @@ -0,0 +1,203 @@ +classdef Beams < handle & matlab.mixin.Copyable + + properties (Access = private) + + BlueBeamDefault = struct('Alias', 'Blue', ... + 'Power', 400e-3, ... + 'Detuning', -1.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 diff --git a/MOT Simulator/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m b/MOT Simulator/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m new file mode 100644 index 0000000..911d39b --- /dev/null +++ b/MOT Simulator/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m @@ -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 diff --git a/MOT Simulator/+Simulator/@Oven/Oven.m b/MOT Simulator/+Simulator/@Oven/Oven.m new file mode 100644 index 0000000..2fa54b2 --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/Oven.m @@ -0,0 +1,177 @@ +classdef Oven < Simulator.MOTCaptureProcess & matlab.mixin.Copyable + + properties (Access = private) + + OvenDefaults = struct('NozzleLength', 60e-3, ... + 'NozzleRadius', 2.60e-3, ... + 'OvenTemperature', 1000); + end + + properties (Access = public) + NozzleLength; + NozzleRadius; + OvenTemperature; + + VelocityCutoff; + ClausingFactor; + ReducedClausingFactor; + ReducedFlux; + end + + properties (Dependent) + ExitDivergence; + Beta; + OvenDistance; + OvenTemperatureinKelvin; + AverageVelocity; + AtomicBeamDensity; + MeanFreePath; + CollisionTime; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods + function this = Oven(varargin) + this@Simulator.MOTCaptureProcess('SimulationMode', '2D', varargin{:}); + this.NozzleLength = this.OvenDefaults.NozzleLength; + this.NozzleRadius = this.OvenDefaults.NozzleRadius; + this.OvenTemperature = this.OvenDefaults.OvenTemperature; + this.ClausingFactor = this.calculateClausingFactor(); + [this.ReducedClausingFactor, ~] = this.calculateReducedClausingFactor(); + end + + function restoreDefaults(this) + this.NozzleLength = this.OvenDefaults.NozzleLength; + this.NozzleRadius = this.OvenDefaults.NozzleRadius; + this.OvenTemperature = this.OvenDefaults.OvenTemperature; + this.ClausingFactor = this.calculateClausingFactor(); + [this.ReducedClausingFactor, ~] = this.calculateReducedClausingFactor(); + end + + end % - lifecycle + + methods + function set.NozzleLength(this,val) + this.NozzleLength = val; + end + function ret = get.NozzleLength(this) + ret = this.NozzleLength; + end + function set.NozzleRadius(this,val) + this.NozzleRadius = val; + end + function ret = get.NozzleRadius(this) + ret = this.NozzleRadius; + end + function set.OvenTemperature(this,val) + this.OvenTemperature = val; + end + function ret = get.OvenTemperature(this) + ret = this.OvenTemperature; + end + function set.VelocityCutoff(this,val) + this.VelocityCutoff = val; + end + function ret = get.VelocityCutoff(this) + ret = this.VelocityCutoff; + end + function set.ClausingFactor(this,val) + this.ClausingFactor = val; + end + function ret = get.ClausingFactor(this) + ret = this.ClausingFactor; + end + function set.ReducedClausingFactor(this,val) + this.ReducedClausingFactor = val; + end + function ret = get.ReducedClausingFactor(this) + ret = this.ReducedClausingFactor; + end + function set.ReducedFlux(this,val) + this.ReducedFlux = val; + end + function ret = get.ReducedFlux(this) + ret = this.ReducedFlux; + end + end % - setters and getters + + methods + function ret = get.Beta(this) + ret = 2 * (this.NozzleRadius/this.NozzleLength); + end + + function ret = get.OvenDistance(this) + ApertureCut = max(2.5e-3,this.NozzleRadius); + ret = (25+12.5)*1e-3 + (this.NozzleRadius + ApertureCut) / tan(15/360 * 2 * pi); + end + + function ret = get.ExitDivergence(this) + Theta_Nozzle = atan((this.NozzleRadius + 0.035/sqrt(2))/this.OvenDistance); % The angle of capture region towards the oven nozzle + Theta_Aperture = 15/360*2*pi; % The limitation angle of the second aperture in the oven + ret = min(Theta_Nozzle,Theta_Aperture); + end + + function ret = get.OvenTemperatureinKelvin(this) + ret = this.OvenTemperature + Helper.PhysicsConstants.ZeroKelvin; + end + + function ret = get.AverageVelocity(this) + %See Background collision probability section in Barbiero + ret = sqrt((8 * pi * Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass)); + end + + function ret = get.AtomicBeamDensity(this) + %See Background collision probability section in Barbiero + ret = this.calculateFreeMolecularRegimeFlux / (this.AverageVelocity * pi * (this.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 diff --git a/MOT Simulator/+Simulator/@Oven/angularDistributionFunction.m b/MOT Simulator/+Simulator/@Oven/angularDistributionFunction.m new file mode 100644 index 0000000..acd3eab --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/angularDistributionFunction.m @@ -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 + diff --git a/MOT Simulator/+Simulator/@Oven/calculateClausingFactor.m b/MOT Simulator/+Simulator/@Oven/calculateClausingFactor.m new file mode 100644 index 0000000..84c7ebc --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/calculateClausingFactor.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m b/MOT Simulator/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m new file mode 100644 index 0000000..92b09dd --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@Oven/calculateReducedClausingFactor.m b/MOT Simulator/+Simulator/@Oven/calculateReducedClausingFactor.m new file mode 100644 index 0000000..919dc91 --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/calculateReducedClausingFactor.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@Oven/initialPositionSampling.m b/MOT Simulator/+Simulator/@Oven/initialPositionSampling.m new file mode 100644 index 0000000..85af394 --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/initialPositionSampling.m @@ -0,0 +1,7 @@ +function ret = initialPositionSampling(this) + n = this.NumberOfAtoms; + phi = 2 * pi * rand(n,1); + rho = this.Beta * 0.5 * this.NozzleLength * sqrt(rand(n,1)); + ret = [-this.OvenDistance * ones(n,1), rho.*cos(phi), rho.*sin(phi)]; +end + diff --git a/MOT Simulator/+Simulator/@Oven/initialVelocitySampling.m b/MOT Simulator/+Simulator/@Oven/initialVelocitySampling.m new file mode 100644 index 0000000..b82775a --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/initialVelocitySampling.m @@ -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 + diff --git a/MOT Simulator/+Simulator/@Oven/velocityDistributionFunction.m b/MOT Simulator/+Simulator/@Oven/velocityDistributionFunction.m new file mode 100644 index 0000000..2f5dc98 --- /dev/null +++ b/MOT Simulator/+Simulator/@Oven/velocityDistributionFunction.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m new file mode 100644 index 0000000..bfdf096 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m @@ -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 diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m new file mode 100644 index 0000000..44d6485 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m @@ -0,0 +1,35 @@ +function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector) + % is the distance between the chamber center and the cross point of push beam and z-axis (along the gravity) + WaveVectorEndPoint = [0, 1, this.DistanceBetweenPushBeamAnd3DMOTCenter/this.PushBeamDistance]; + WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint); + Origin=[0,0,0]; + + PushBeamObj = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), this.Beams)}; + PushBeamDetuning = PushBeamObj.Detuning; + PushBeamRadius = PushBeamObj.Radius; + PushBeamWaist = PushBeamObj.Waist; + PushBeamWaveNumber = PushBeamObj.WaveNumber; + PushBeamLinewidth = PushBeamObj.Linewidth; + PushBeamSaturationParameter = 0.25 * PushBeamObj.SaturationParameter; + + SaturationIntensity = this.calculateLocalSaturationIntensity(PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, PushBeamRadius, PushBeamWaist); + + DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * PushBeamWaveNumber; + + Detuning = PushBeamDetuning - DopplerShift; + + s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./PushBeamLinewidth).^2)); + a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(PushBeamLinewidth * 0.5); + a_push = a_sat .* (s_push/(1+s_push)); + + if this.SpontaneousEmission + a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, PushBeamLinewidth, PushBeamWaveNumber); + else + a_scatter = [0,0,0]; + end + + a_total = a_push + a_scatter; + + ret = a_total(1:3); +end + diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m new file mode 100644 index 0000000..9cf4a31 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m @@ -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 + diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m new file mode 100644 index 0000000..a06802a --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m new file mode 100644 index 0000000..074fad7 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m @@ -0,0 +1,22 @@ +function ret = calculateCaptureVelocity(this, ovenObj, PositionVector, VelocityVector) + VelocityUnitVector = VelocityVector./norm(VelocityVector); + UpperLimit = 500; + LowerLimit = 0; + + for Index = 1:500 + InitialVelocity = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; + ParticleDynamicalQuantities = this.solver(PositionVector, InitialVelocity); + FinalPositionVector = ParticleDynamicalQuantities(end, 1:3); + if rssq(FinalPositionVector) <= ovenObj.OvenDistance + LowerLimit = 0.5 * (UpperLimit + LowerLimit); + else + UpperLimit = 0.5 * (UpperLimit + LowerLimit); + end + + if UpperLimit - LowerLimit < 1 + ret = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector; + break; + end + end +end + diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m new file mode 100644 index 0000000..eada4cb --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m @@ -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 diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m new file mode 100644 index 0000000..5bc099a --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m new file mode 100644 index 0000000..68513eb --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m new file mode 100644 index 0000000..39b0014 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m new file mode 100644 index 0000000..6e1aef0 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/exitCondition.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/exitCondition.m new file mode 100644 index 0000000..33c1406 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/exitCondition.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m new file mode 100644 index 0000000..abca489 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m new file mode 100644 index 0000000..4fc4e1c --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m @@ -0,0 +1,8 @@ +function ret = magneticFieldForMOT(this, r) + ret = zeros(1,4); + alpha = this.MagneticGradient; + ret(1) = r(3)*alpha; + ret(2) = 0; + ret(3) = r(1)*alpha; + ret(4) = sqrt(ret(1)^2+ret(2)^2+ret(3)^2); +end \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/runSimulation.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/runSimulation.m new file mode 100644 index 0000000..f26b3b1 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/runSimulation.m @@ -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 \ No newline at end of file diff --git a/MOT Simulator/+Simulator/@TwoDimensionalMOT/solver.m b/MOT Simulator/+Simulator/@TwoDimensionalMOT/solver.m new file mode 100644 index 0000000..a2c3cf1 --- /dev/null +++ b/MOT Simulator/+Simulator/@TwoDimensionalMOT/solver.m @@ -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 + diff --git a/MOT Simulator/test_MOTCaptureProcessSimulation.m b/MOT Simulator/test_MOTCaptureProcessSimulation.m new file mode 100644 index 0000000..2c027ad --- /dev/null +++ b/MOT Simulator/test_MOTCaptureProcessSimulation.m @@ -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(); +