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