% This function displays some parameters after tracking
function disp_results(beam)
    format long;
    emitt = placet_get_emittance(beam);
    ex = emitt(1)
    ey = emitt(2)
    disp("");
    twiss = placet_get_twiss_matrix(beam);
    beta_x  =  twiss(1,1)
    beta_y  =  twiss(3,3)
    alpha_x = -twiss(1,2)
    alpha_y = -twiss(3,4)
    disp("");
    x  = mean(beam(:,2))
    y  = mean(beam(:,3))
    xp = mean(beam(:,5))
    yp = mean(beam(:,6))
    disp("");
    z  = mean(beam(:,4))
    sz =  std(beam(:,4))
    disp("");
    E  = mean(beam(:,1))
    dE =  std(beam(:,1))
    disp("");
endfunction

% This funciton returns the beam and Twiss parameters at the end of a section
function results = get_tracking_results(beam, sigma_cut)

    % Cut the particles outside 3 sigma before calculating beam parameters
    beam = cut_distribution(beam, sigma_cut);

    % Beam parameters
    emittance = placet_get_emittance(beam);
    ex  = emittance(1);
    ey  = emittance(2);
    x   = mean(beam(:,2));
    y   = mean(beam(:,3));
    z   = mean(beam(:,4));
    xp  = mean(beam(:,5));
    yp  = mean(beam(:,6));
    sx  = std(beam(:,2));
    sy  = std(beam(:,3));
    sz  = std(beam(:,4));
    sxp = std(beam(:,5));
    syp = std(beam(:,6));
    E   = mean(beam(:,1));
    dE  = std(beam(:,1));

    % Twiss parameters
    twiss_matrix = placet_get_twiss_matrix(beam);
    bx =  twiss_matrix(1,1);
    by =  twiss_matrix(3,3);
    ax = -twiss_matrix(1,2);
    ay = -twiss_matrix(3,4);

    # Collect the parameters into a single matrix
    results = [ex ey x y z xp yp sx sy sz sxp syp E dE bx by ax ay];

endfunction    

% This function writes the results of each section to a file
function save_tracking_results(filename, rtml, ml, bds)
    results = [rtml; ml; bds];
    save("-text", filename, "results");
endfunction
