%% 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();