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

    for i=1:length(Y)
        plot(DynamicalQuantities(i,:,1),DynamicalQuantities(i,:,4),'w','linewidth',1.3)
    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