Added code to allow plotting of Local saturation Intensity distribution for the beams.

test_MOTCaptureProcessSimulation.m

@@ -75,15 +75,16 @@ Plotter.plotAngularDistributionForDifferentBeta(Oven, Beta)
 Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot!
 
 %% - Plot Phase Space with Acceleration Field
-MOT2D.SidebandBeam        = true;
+MOT2D.SidebandBeam        = false;
+MOT2D.MagneticGradient    = 0.4;
 CoolingBeam               = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
-CoolingBeam.Power         = 0.2;
+CoolingBeam.Power         = 0.3;
-CoolingBeam.Detuning      = -1.5*Helper.PhysicsConstants.BlueLinewidth;
+CoolingBeam.Detuning      = -1.64*Helper.PhysicsConstants.BlueLinewidth;
-CoolingBeam.Waist         = 16e-03;
+CoolingBeam.Waist         = 15e-03;
 SidebandBeam              = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
-SidebandBeam.Power        = 0.6;
+SidebandBeam.Power        = 0.5;
-SidebandBeam.Detuning     = -5*Helper.PhysicsConstants.BlueLinewidth;
+SidebandBeam.Detuning     = -4*Helper.PhysicsConstants.BlueLinewidth;
-SidebandBeam.Waist        = 16e-03;
+SidebandBeam.Waist        = 15e-03;
 MOT2D.NumberOfAtoms       = 50;
 MinimumVelocity           = 0;
 MaximumVelocity           = 150;
@@ -265,8 +266,9 @@ clear OptionsStruct
 MOT2D.NumberOfAtoms           = 10000;
 MOT2D.SidebandBeam            = false;
 MOT2D.PushBeam                = false;
+MOT2D.BackgroundCollision     = false;
 CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
-CoolingBeam.Power             = 0.7;
+CoolingBeam.Power             = 0.4;
 NumberOfPointsForFirstParam   = 10; %iterations of the simulation
 NumberOfPointsForSecondParam  = 10;
 NumberOfPointsForThirdParam   = 6;
@@ -300,3 +302,117 @@ options                                             = Helper.convertstruct2cell(
 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 = 1000;
+t = 2*pi*rand(n,1);
+r = BeamRadius*sqrt(rand(n,1));
+% xmin = -15e-03;
+% xmax = 15e-03;
+% x    = xmin+rand(n,1)*(xmax-xmin);
+x = r.*cos(t);
+y = ones(n,1) * 2e-03;
+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 = 100;
+columns = 100;
+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(-BeamWaist, BeamWaist, n), linspace(-BeamWaist, BeamWaist, 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();
+