Calculations/+Plotter/plotInitialVeloctiySamplingVsAngle.m

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