/**
    Exosphere simulation code as described in Kegerreis et al. (2016) and
    the README.pdf documentation.

    Publically available at icc.dur.ac.uk/index.php?content=Research/Topics/O13

    Simulation function.
**/


/*
    0. Includes etc.
    1. Constants
    2. Simulation Function
*/


////////////////////////////////////////////////////////////////////////////////
//                              0. Includes etc.                              //
////////////////////////////////////////////////////////////////////////////////


#include "functions.h"
#include "simulation.h"
#include "nr_cpp/nr.h"

#include <iostream>
#include <cmath>
#include <numeric>
#include <ctime>
#include <sys/stat.h>
#include <string>

#include <mpi.h>

// to_string() not working so
#include <sstream>
#define SSTR( x ) dynamic_cast< std::ostringstream & >( \
        ( std::ostringstream() << std::dec << x ) ).str()

using namespace std;


////////////////////////////////////////////////////////////////////////////////
//                              1. Constants                                  //
////////////////////////////////////////////////////////////////////////////////


// Standard
const double pi     = 3.14159265358979323846;
const double conv   = pi/180;           // Degrees to radians (rad/deg)
const double k      = 1.3806485e-23;    // Boltzmann constant (J K^-1)
const double G      = 6.67384e-11;      // Gravitational const. (m^3 kg^-1 s^-2)
const double h      = 6.626176e-34;     // Planck's constant (J s)
const double amu    = 1.66053892e-27;   // Atomic mass unit (kg)
const double year   = 365.25*24*60*60;  // Earth year (s) (= 12.37 lunar days)
// Moon
const double R      = 1738e3;           // Radius of the Moon (m)
const double M      = 7.34767e22;       // Mass of the Moon (kg)
const double day    = 2551443;          // Lunar Synodic day (s)
const double g_0    = G*M/(R*R);        // Surface gravity = 1.6234 (m s^-2)
const double v_esc  = sqrt(2*G*M/R);    // Escape velocity = 2375.49 (m s^-1)


////////////////////////////////////////////////////////////////////////////////
//                              2. Simulation Function                        //
////////////////////////////////////////////////////////////////////////////////


void run_simulation(const int num_run, const int num_particle,
                    const int percent, const double dt, const double start_time,
                    const double end_time, const int start_N_bounce,
                    const int end_N_bounce, const double mass, const double tau,
                    const double m, double Q, double P_adsorb,
                    double P_squirrel, double f_N, double f_S,
                    string source, string emerge_dist, string losses,
                    string T_dist, string cold_traps, string record_dist,
                    string surf_intn, string misc, string rates,
                    string backup_results,
                    int job_ID, int job_index)
{
    /// Performs the simulation as described in Kegerreis et al. (2016) and \
        the README.pdf documentation.

    // Check the string arguments from main.cpp, print error if invalid
    check_arguments(source, emerge_dist, losses, T_dist, cold_traps,
                    record_dist, surf_intn, misc, rates, backup_results);

    //=========================
    // Default Simulation Variables
    //=========================


    //========
    // Basic Parameters
    //========
    double  theta_0, phi_0, theta, phi, // Latitude and longitude (rad)
            theta_source, phi_source,   // Source coordinates (rad)
            phi_,                       // Local time of day (rad)
            time,                       // Simulation time (s)
            T, T_min, T_0, T_1,         // Temperatures (K)
            v_0,                        // Initial speed (m s^-1)
            t,                          // Time of flight (s)
            // Hop geometry parameters as defined in the paper, appendix A
            alpha, eta, h_max, a, e, delta, epsilon,
            t_res,                      // Residence time (s)
            flux,                       // Particle flux (contribution to the \
                                           density) ((m s^-1)^-1)
            N_bounce,                   // Number of bounces
            phi_offset,                 // Longitude offset for the T map (rad)
            xi_1, xi_2, xi_3, xi_4,     // Random numbers
            Prob,                       // Probability (mutli-use)
            sig_offset,                 // Longitude offsets sigma (rad)
            // Fractions of particles lost by each mechanism
            frac_captured, frac_destroyed, frac_escaped,
            // Means and standard deviations of ^
            mean_frac_captured, mean_frac_destroyed, mean_frac_escaped,
            stdev_captured, stdev_destroyed, stdev_escaped;

    int     // Number of particles lost by each mechanism
            N_captured, N_destroyed, N_escaped, N_ended,
            // Indices for accessing a specific array element
            theta_index, phi__index, phi_index;


    //========
    // Misc
    //========
    // Date and time
    const string    date_time   = get_date_time('y', 'y');
    clock_t         t_start, t_end;

    // Main path
    const string    dir         = "/cosma/home/jkeger/exosphere/";
    // Data files path
    const string    dir_data    = "/cosma/home/jkeger/exosphere/data/";

    // Simulation particle source rate (s^-1)
    const double    N_sim   = num_run*num_particle / end_time;


    //========
    // Arrays
    //========
    // phi/phi_ and theta array dimensions
    int num_phi     = 360,
        num_phi_    = num_phi,
        num_theta   = num_phi/2;
    // phi/phi_ and theta maximum values (rad)
    double  max_phi     = 2 * pi,
            max_phi_    = 2 * pi,
            max_theta   = pi;

    // Store loss results from each simulation run
    vector<double>  A1_frac_captured, A1_frac_destroyed, A1_frac_escaped;

    // Resize the arrays
    A1_frac_captured.resize(num_run);
    A1_frac_destroyed.resize(num_run);
    A1_frac_escaped.resize(num_run);

    vector< vector<double> >    // Velocity generation array
                                A2_Maxwell_cdf,
                                // Temperature map (generated for each run)
                                A2_H_temperature,
                                // Temperature map longitude offsets
                                A2_longitude_offset;

    // Set temperature parameters
    if ((T_dist == "H_T_simple") || (T_dist == "H_T_offsets")) {
        T_min   = 60;
        T_0     = 130;
        T_1     = 262;
    }
    if (T_dist == "H_T_offsets") {
        sig_offset  = 4.5 * conv;
    }

    // Create array of Maxwell cdfs
    if (emerge_dist == "M_B") {
        A2_Maxwell_cdf  = make_A2_Maxwell_cdf(T_0, T_1, T_min, m);
    }
    else if (emerge_dist == "M_B_flux") {
        A2_Maxwell_cdf  = make_A2_Maxwell_flux_cdf(T_0, T_1, T_min, m);
    }


    //=========================
    // Variables for some options in main.cpp
    //=========================


    //========
    // source = "local_source" or "cont_local_source"
    //========
    double  // Source latitude and longitude boundary coordinates (deg)
            theta_l_1, theta_l_2, phi_l_1, phi_l_2,
            // Boundary random numbers that would give the boundary coordinates
            xi_phi_1, xi_phi_2, xi_theta_1, xi_theta_2;

    if ((source == "local_source") || (source == "cont_local_source")) {
        // In degrees for nicer printing. Only the xis are actually used.
        theta_l_1   = 89;
        theta_l_2   = 90;
        phi_l_1     = 339;
        phi_l_2     = 340;

        xi_phi_1    = 0.5 * phi_l_1*conv / pi;
        xi_phi_2    = 0.5 * phi_l_2*conv / pi;
        xi_theta_1  = 0.5*(1 - cos(theta_l_1*conv));
        xi_theta_2  = 0.5*(1 - cos(theta_l_2*conv));
    }


    //========
    // source = "K_source" or "eqm_source"
    //========
    int K_index, eqm_index;

    // Index for finding the right element from a random number.
    vector<double>::iterator    bound_index;

    // Load in data for "K_source" LPGRS potassium source map cdf
    string  data0 = dir_data, data1 = dir_data, data2 = dir_data,
            data3 = dir_data, data4 = dir_data;

    ifstream                    file0(data0.append("K_ppm.txt"));
    // Instead, use a top-hat source from a 2000 ppm cut on the LPGRS map
//    ifstream                    file0(data0.append("K_ppm_cut.txt"));
    istream_iterator<double>    start0(file0), end;
    vector<double>              A1_K_ppm(start0, end);

    ifstream                    file1(data1.append("K_theta_0.txt"));
    istream_iterator<double>    start1(file1);
    vector<double>              A1_K_theta_0(start1, end);

    ifstream                    file2(data2.append("K_theta_1.txt"));
    istream_iterator<double>    start2(file2);
    vector<double>              A1_K_theta_1(start2, end);

    ifstream                    file3(data3.append("K_phi_0.txt"));
    istream_iterator<double>    start3(file3);
    vector<double>              A1_K_phi_0(start3, end);

    ifstream                    file4(data4.append("K_phi_1.txt"));
    istream_iterator<double>    start4(file4);
    vector<double>              A1_K_phi_1(start4, end);

    // Load in equilibrium source cdf for "eqm_source"
    string                      dir_eqm = dir_data;
    ifstream                    file_eqm(dir_eqm.append("eqm_source.txt"));
    istream_iterator<double>    start_eqm(file_eqm);
    vector<double>              A1_eqm_source(start_eqm, end);


    //========
    // source = "transient_source"
    //========
    // Location and time of the transient source (rad), (s)
    double  theta_trans, phi_trans, time_trans;

    theta_trans = 70 * conv;
    phi_trans   = 150 * conv;
    time_trans  = 270/360 * day;


    //========
    // record_dist = "density" or "density_stuck"
    //========
    // Array to record the density distribution of particles by latitude and \
       local time of day
    vector< vector<double> >    A2_density;

    // Resize the array if it will be used
    if ((record_dist == "density") || (record_dist == "density_stuck")) {
        A2_density.resize(num_theta);
        for (int i_theta=0; i_theta<A2_density.size(); i_theta++) {
            A2_density[i_theta].resize(num_phi_);
        }
    }

    // For record_dist = "density_stuck"
    double  t_stuck,                    // Time residing on the surface (s)
            phi__end;                   // phi_ at the simulation's end (rad)
    int     phi__index_stuck;


    //========
    // record_dist = "densities" or "all_densities"
    //========
    int num_time, num_day, time_index;

    // Array to record the density distribution of particles in one of two ways
    vector< vector< vector<double> > >  A3_density;

    // Resize the array if it will be used
    // One map of latitude and longitude for each local time of day
    if (record_dist == "densities") {
        A3_density.resize(num_phi_);
        // Resize each empty map
        for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
            A3_density[i_phi_].resize(num_theta);
            for (int i_theta=0; i_theta<num_theta; i_theta++) {
                A3_density[i_phi_][i_theta].resize(num_phi);
            }
        }
    }
    // One map of latitude and local time of day for each step in time
    else if (record_dist == "all_densities") {
        // Number of maps to make (num_time for each day)
        num_time    = 360;
        num_day     = floor((end_time-start_time)/day) + 1;
        // One map at each time with timestep of day/num_time
        A3_density.resize(num_time * num_day);
        // Resize each empty map
        for (int i_day=0; i_day<num_day; i_day++) {
            for (int i_t=0; i_t<num_time; i_t++) {
                A3_density[(i_day*num_time) + i_t].resize(num_theta);
                for (int i_theta=0; i_theta<num_theta; i_theta++) {
                    A3_density[(i_day*num_time)+i_t][i_theta].resize(num_phi_);
                }
            }
        }
    }


    //========
    // record_dist = "at_altitude"
    //========
    int alt,                            // Altitude (m)
        alt_step,                       // Distance between each altitude (m)
        alt_step_max,                   // Maximum altitude (m)
        num_alt;

    double  v,                          // Speed at altitdue (m s^-1)
            // Extra coordinates for the positions at altitude (rad)
            theta_1, theta_2, phi_1, phi_2, phi__1, phi__2,
            delta_p,                    // delta angle to the peak at alt. (rad)
            delta_1, delta_2,           // deltas to the positions at alt. (rad)
            t_p;                        // Time of flight to the peak (s)


    // Array to record all particles at altitude so that the mean and std can \
       be taken afterwards
    vector< vector< vector< vector<double> > > >    A4_density_alt;

    // Arrays to hold the mean and variance of the density at altitude
    vector< vector< vector<double> > >  A3_density_alt, A3_variance_alt;

    // Resize the arrays if they will be used
    if (record_dist == "at_altitude") {
        // Normal theta, phi_, then altitude third dimension
        num_alt     = 15;
        alt_step    = 10 * 1e3;
        A4_density_alt.resize(num_theta);
        A3_density_alt.resize(num_theta);
        A3_variance_alt.resize(num_theta);

        for (int i_theta=0; i_theta<num_theta; i_theta++) {
            A4_density_alt[i_theta].resize(num_phi_);
            A3_density_alt[i_theta].resize(num_phi_);
            A3_variance_alt[i_theta].resize(num_phi_);

            for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                A4_density_alt[i_theta][i_phi_].resize(num_alt);
                A3_density_alt[i_theta][i_phi_].resize(num_alt);
                A3_variance_alt[i_theta][i_phi_].resize(num_alt);
            }
        }
    }


    //========
    // cold_traps = "local_traps"
    //========
    // Boundary coordinates of the north and south traps (rad)
    double  theta_trap_N_1  = 0 * conv,
            theta_trap_N_2  = 15 * conv,
            phi_trap_N_1    = 0 * conv,
            phi_trap_N_2    = 10 * conv,

            theta_trap_S_1  = 165 * conv,
            theta_trap_S_2  = 180 * conv,
            phi_trap_S_1    = 0 * conv,
            phi_trap_S_2    = 10 * conv;


    //========
    // cold_traps = "seasonal"
    //========
    double  f_season_N  = //8e-3,         // Peak seasonal trap fractional area
                          //4e-2,
                          1.8e-1,
            f_season_S  = f_season_N;


    //========
    // surf_intn = "squirrelling"
    //========
    double  t_sq,                       // Time spent squirrelling (s)
            // Boundary phi_ for squirrelling in and out of the regolith (rad)
            phi__sq_in_1    = 90 * conv,
            phi__sq_in_2    = 270 * conv,
            phi__sq_out_1   = 270 * conv,
            // Add 360 if phi__sq_out_2 is past noon so that (_out - _in) > 0
            phi__sq_out_2   = (90 + 360) * conv;


    //========
    // misc = "local_Q" or "mare_T_offset"
    //========
    double  // Mare region boundary coordinates (rad)
            theta_m_1   = 30 * conv,
            theta_m_2   = 110 * conv,
            phi_m_1     = 300 * conv,
            phi_m_2     = 350 * conv;

    // For misc = "local_Q"
    double  Q_m = 26,                   // Mare desorption energy (kJ mol^-1)
            Q_h = 28;                   // Highland desorp. energy (kJ mol^-1)

    // For misc = "mare_T_offset"
    double  T_mare  = 10;               // Temperature offset for the mare (K)


    //========
    // misc = "solar_flare"
    //========
    double  t_flare_1   = 200*day,      // Start and end times of the flare (s)
            t_flare_2   = 200*day + 60*60;


    //========
    // rates = "yes_rates" or "LADEE_lat_rates"
    //========
    int num_time_rates  = 1e4;

    vector<double>  // Simulation times (s) of source and loss for each particle
                    A1_time_start, A1_time_capt, A1_time_dest,
                    A1_time_esc,
                    // and for the seasonal cold traps
                    A1_time_seasonal_capt_N, A1_time_seasonal_rel_N,
                    A1_time_seasonal_capt_S, A1_time_seasonal_rel_S;

    vector<int>     // Binned number of particles by their source and loss times
                    A1_bin_time_start, A1_bin_time_capt, A1_bin_time_dest,
                    A1_bin_time_esc,
                    // and for the seasonal cold traps
                    A1_bin_time_seasonal_capt_N, A1_bin_time_seasonal_rel_N,
                    A1_bin_time_seasonal_capt_S, A1_bin_time_seasonal_rel_S;


    // ================================================================== //


    // Make something dependent on job_index for easy multi-run jobs
    bool    use_job_index   = true;

    // Set different parameter values for each array job
    if (use_job_index) {
        if (job_index != 0) {
//            Q   = 23 + job_index;
//            P_squirrel  = 1e-3 * (2*job_index - 1);

            f_N = pow(10, -1.25 - 0.25*job_index);
            f_S = f_N;

//            f_season_N  = job_index * 1e-2 + 1e-1;
//            f_season_S  = f_season_N;
        }
    }


    // ================================================================== //
    MPI_Init(NULL, NULL);
    // ======================== Start MPI


    //========
    // Get total number of processors and rank of this processor
    //========
    int world_size, world_rank;
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);
    MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

    // Get number of runs for each processor, print error if not integer
    int runs_per_proc   = num_run / world_size;
    if (num_run%world_size != 0) {
        cout << "Warning: number of runs not divisible by number of processors!"
             << endl;
    }

    // Resize vectors which depend on runs_per_proc, if used
    if (record_dist == "at_altitude") {
        for (int i_theta=0; i_theta<num_theta; i_theta++) {
            for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                for (int i_alt=0; i_alt<num_alt; i_alt++) {
                    A4_density_alt[i_theta][i_phi_][i_alt].resize(
                            runs_per_proc);
                }
            }
        }
    }


    //========
    // Set up RNG for each processor
    //========
    // Generate a unique seed for each processor's RNG
    int seed, A1_seed[world_size];


    // Use the job ID (include the job index)
    int idum    = - job_ID - job_index;
    // Or a manual seed instead
//    int idum    = - 123456;


    // Print initial seed
    if (world_rank == 0) {
        cout << "RNG seed = " << idum << endl;
    }

    // Make a list of the seeds for all processors
    for (int i_rank=0; i_rank<world_size; i_rank++) {
        A1_seed[i_rank] = NR::ran3(idum) * -1e6;
    }

    // Set each rank's seed
    seed    = A1_seed[world_rank];


    //========
    // [0]: Initial printing
    //========
    // Print all the relevant input parameters and information
    if (world_rank == 0) {
        cout << endl << "//=========================" << endl
             << "    Exosphere simulation " << endl
             << "    as described in Kegerreis et al. (2016) " << endl
             << "    and the README.pdf documentation." << endl
             << "//=========================" << endl << endl
             << "Running on " << world_size << " processors." << endl
             << "Job ID <" << job_ID << ">";
        if (job_index != 0) {
            cout << "[" << job_index << "]";
        }
        cout <<  endl << endl;

        // Print the input parameters
        print_debug3(num_run, runs_per_proc, num_particle);
        cout << endl;
        print_debug6(percent, mass, tau, N_sim, Q, P_adsorb);
        print_debug2(f_N, f_S);
        cout << endl;
        // Start and end times in units of lunar days for nicer printing
        double  end_day = end_time/day, start_day = start_time/day;
        print_debug6(start_N_bounce, end_N_bounce, start_day, end_day, num_phi,
                     num_theta);
        cout << endl;
        print_debug6(source, emerge_dist, losses, T_dist, cold_traps,
                     record_dist);
        print_debug3(surf_intn, rates, backup_results);
        cout << endl;
        if (source == "local_source") {
            print_debug4(theta_l_1, theta_l_2, phi_l_1, phi_l_2);
        }
        if (cold_traps == "seasonal") {
            print_debug2(f_season_N, f_season_S);
        }
        if (record_dist == "all_densities") {
            print_debug(num_time);
        }
        if (surf_intn == "squirrelling") {
            print_debug(P_squirrel);
        }
        if (misc == "local_Q") {
            print_debug2(Q_m, Q_h);
        }
        else if (misc == "mare_T_offset") {
            print_debug(T_mare);
        }
        else if (misc == "solar_flare") {
            print_debug2(t_flare_1, t_flare_2);
        }
        cout << endl;

        if (use_job_index) {
            cout << "Using job_index to set parameters!" << endl << endl;
        }
        else {
            cout << "Not using job_index to set parameters." << endl << endl;
        }


        // Print the start time
        cout << "Started at ";
        print_date_time('y', 'y');
        cout << endl;
    }


    // Set start time just before the main simulation loops start
    t_start = clock();


    // =================================
    // Simulation Runs
    // =================================
    for (int i_run=0; i_run<runs_per_proc; i_run++) {
        // Reset loss counters for each run
        N_captured  = 0;
        N_destroyed = 0;
        N_escaped   = 0;
        N_ended     = 0;

        // Print progress every few runs for every processor
        print_progress_MPI(i_run, runs_per_proc, percent, t_start, world_rank,
                           world_size);

        // Set temperature and longitude offset maps
        if (T_dist == "H_T_simple") {
            A2_H_temperature    = make_A2_H_temperature(T_0, T_1, num_phi_,
                                                         num_theta);
        }
        else if (T_dist == "H_T_offsets") {
            A2_longitude_offset = make_A2_longitude_offset(seed - 1, sig_offset,
                                                           num_phi_, num_theta);
            A2_H_temperature    = make_A2_H_temperature(T_0, T_1, num_phi_,
                                                        num_theta);
        }


        // =================================
        // Particle Runs
        // =================================
        for (int i_par=0; i_par<num_particle; i_par++) {
            // Reset simulation time and N_bounce for each particle
            time        = 0;
            N_bounce    = 0;


            //========
            // Generate starting coordinates
            //========
            xi_1    = NR::ran3(seed);
            xi_2    = NR::ran3(seed);

            if (source == "iso_source") {
                theta_source    = rand_theta_iso(xi_1);
                phi_source      = rand_phi_iso(xi_2);
            }
            else if (source == "solar_source") {
                xi_3    =   NR::ran3(seed);

                theta_source    = rand_theta_solar(xi_1, xi_2);
                phi_source      = rand_phi_solar(xi_1, xi_2, xi_3);
            }
            else if (source == "local_source") {
                // Scale xi's to be between xi_theta,phi_1,2 then use iso_source
                xi_1    = xi_1 * (xi_theta_2-xi_theta_1) + xi_theta_1;
                xi_2    = xi_2 * (xi_phi_2-xi_phi_1) + xi_phi_1;

                theta_source    = rand_theta_iso(xi_1);
                phi_source      = rand_phi_iso(xi_2);
            }
            else if (source == "cont_iso_source") {
                theta_source    = rand_theta_iso(xi_1);
                phi_source      = rand_phi_iso(xi_2);

                // Set random start time
                xi_3    = NR::ran3(seed);
                time    = xi_3*end_time;
            }
            else if (source == "cont_local_source") {
                // Scale xi's to be between xi_theta,phi_1,2 then use iso_source
                xi_1    = xi_1 * (xi_theta_2-xi_theta_1) + xi_theta_1;
                xi_2    = xi_2 * (xi_phi_2-xi_phi_1) + xi_phi_1;

                theta_source    = rand_theta_iso(xi_1);
                phi_source      = rand_phi_iso(xi_2);

                // Set random start time
                xi_3    = NR::ran3(seed);
                time    = xi_3*end_time;
            }
            else if (source == "K_source") {
                // One random number needed to select the bin.
                // Find the index of the first value larger than xi
                bound_index = lower_bound(A1_K_ppm.begin(), A1_K_ppm.end(),
                                          xi_1);
                K_index     = (bound_index - A1_K_ppm.begin());

                // Like local_source, get random location within the \
                   corresponding bin.
                // Boundary xi's that would give the boundary bin coords
                xi_phi_1    = 0.5 * A1_K_phi_0[K_index] * conv/pi;
                xi_phi_2    = 0.5 * A1_K_phi_1[K_index] * conv/pi;
                xi_theta_1  = 0.5 * (1 - cos( A1_K_theta_0[K_index] * conv));
                xi_theta_2  = 0.5 * (1 - cos( A1_K_theta_1[K_index] * conv));

                // Scale xi's to be between xi_theta,phi_1,2 then use iso_source
                xi_3    = NR::ran3(seed);
                xi_2    = xi_2 * (xi_theta_2-xi_theta_1) + xi_theta_1;
                xi_3    = xi_3 * (xi_phi_2-xi_phi_1) + xi_phi_1;

                theta_source    = rand_theta_iso(xi_2);
                phi_source      = rand_phi_iso(xi_3);

                // Shift phi_source by -180 degrees because K map has shifted 0
                phi_source  -= pi;
                if (phi_source < 0) {
                    phi_source  += 2*pi;
                }

                // Set random start time
                xi_3    = NR::ran3(seed);
                time    = xi_3*end_time;
            }
            else if (source == "transient_source") {
                // One-time point source
                theta_source    = theta_trans;
                phi_source      = phi_trans;
                time            = time_trans;
            }
            else if (source == "eqm_source") {
                // One random number needed to select the bin.
                // Find the index of the first value larger than xi
                bound_index = lower_bound(A1_eqm_source.begin(),
                                          A1_eqm_source.end(), xi_1);
                eqm_index   = (bound_index - A1_eqm_source.begin());

                // Get theta and phi from the random index
                theta_1 = floor(eqm_index / 360);
                theta_2 = theta_1 + 1;
                phi_1   = eqm_index - theta_1 * 360;
                phi_2   = phi_1 + 1;

                // Like local source, get random location within the \
                   corresponding bin.
                // Boundary xi's that would give the boundary bin coords
                xi_phi_1    = 0.5 * phi_1 * conv/pi;
                xi_phi_2    = 0.5 * phi_2 * conv/pi;
                xi_theta_1  = 0.5 * (1 - cos(theta_1 * conv));
                xi_theta_2  = 0.5 * (1 - cos(theta_2 * conv));

                // Scale xi's to be between xi_theta,phi_1,2 then use iso_source
                xi_3    = NR::ran3(seed);
                xi_2    = xi_2 * (xi_theta_2-xi_theta_1) + xi_theta_1;
                xi_3    = xi_3 * (xi_phi_2-xi_phi_1) + xi_phi_1;

                theta_source    = rand_theta_iso(xi_2);
                phi_source      = rand_phi_iso(xi_3);
            }

            // Initialise coordinates for the first hop
            theta_0 = theta_source;
            phi_0   = phi_source;
            theta   = theta_source;
            phi     = phi_source;

            // Record when the particle enters the simulation
            if (rates == "yes_rates") {
                A1_time_start.push_back(time);
            }
            else if (rates == "LADEE_lat_rates") {
                if (is_LADEE_lat(theta)) {
                    A1_time_start.push_back(time);
                }
            }


            // =================================
            // Migrate until loss
            // =================================
            while (true) {
                // Count number of bounces
                N_bounce++;


                //========
                // Surface Interactions
                //========
                // Probability of adsorption
                xi_1    = NR::ran3(seed);
                // Skip surface interactions if the particle does not adsorb
                if (xi_1 < P_adsorb) {
                    // Squirrelling down to reappear later
                    if (surf_intn == "squirrelling") {
                        if ((phi__sq_in_1 < phi_) && (phi_ < phi__sq_in_2)) {
                            xi_1    = NR::ran3(seed);
                            xi_2    = NR::ran3(seed);
                            // If squirrelled
                            if (xi_1 < P_squirrel) {
                                // Advance time for Moon to rotate until reentry
                                t_sq    = ((phi__sq_out_1 + xi_2 *
                                           (phi__sq_out_2 - phi__sq_out_1))
                                           - phi_) * day/(2*pi);
                                time    += t_sq;
                            }
                        }
                    }

                    // Residence time
                    if (surf_intn != "zero") {
                        // Record when the particle first sticks
                        if (record_dist == "density_stuck") {
                            t_stuck             = time;
                            phi__index_stuck    = get_index(phi_, max_phi_,
                                                            num_phi_);
                        }

                        // Check every dt for unsticking until it desorbs
                        while (true) {
                            // Get local temperature
                            if (T_dist == "H_T_simple") {
                                // Find the phi_ and temperature
                                phi_        = phi_prime(phi, time);
                                phi__index  = get_index(phi_, max_phi_,
                                                        num_phi_);
                                theta_index = get_index(theta, max_theta,
                                                        num_theta);

                                T   = A2_H_temperature[phi__index][theta_index];
                            }
                            else if (T_dist == "H_T_offsets") {
                                // First find and use the phi offset
                                phi_index   = get_index(phi, max_phi, num_phi);
                                theta_index = get_index(theta, max_theta,
                                                        num_theta);

                                phi_offset  = A2_longitude_offset[phi_index]
                                                                 [theta_index];

                                // Find phi_ and temperature for that offset phi
                                phi_        = phi_prime(phi + phi_offset, time);
                                phi__index  = get_index(phi_, max_phi_,
                                                        num_phi_);

                                T   = A2_H_temperature[phi__index][theta_index];
                            }

                            // Mare offset temperature
                            if (misc == "mare_T_offset") {
                                if (is_in_mare(theta, phi, theta_m_1, theta_m_2,
                                               phi_m_1, phi_m_2)) {
                                    T   += T_mare;
                                }
                            }

                            // Get probability to unstick
                            if (misc == "local_Q") {
                                // Desorption energy is location dependent
                                if (is_in_mare(theta, phi, theta_m_1, theta_m_2,
                                               phi_m_1, phi_m_2)) {
                                    Q   = Q_m;
                                }
                                else {
                                    Q   = Q_h;
                                }
                            }

                            t_res   = residence_time(T, Q);
                            Prob    = exp(-dt/t_res);
                            xi_1    = NR::ran3(seed);

                            // Break if it unsticks
                            if (xi_1 > Prob) {
                                // Keep residence time != 0 if t_res is small
                                if (t_res < dt) {
                                    time    += t_res;
                                }

                                // Record the time stuck on the surface
                                if (record_dist == "density_stuck") {
                                    if ((N_bounce > start_N_bounce) &&
                                        (time > start_time)) {
                                        // Get indices of relevant array bin
                                        theta_index = get_index(theta,
                                                                max_theta,
                                                                num_theta);

                                        // Record particle's time on the surface
                                        A2_density[theta_index]
                                                  [phi__index_stuck]    +=
                                                time - t_stuck;
                                    }
                                }

                                break;
                            }
                            // Advance time by dt and repeat
                            time    += dt;

                            // If still stuck, check if the phi_ bin has changed
                            if (record_dist == "density_stuck") {
                                phi_        = phi_prime(phi_0, time);
                                phi__index  = get_index(phi_, max_phi_,
                                                        num_phi_);

                                // Record the particle in the bin it just left
                                if (phi__index != phi__index_stuck) {
                                    if (is_recording_time(time, N_bounce,
                                                          start_time,
                                                          start_N_bounce)) {
                                        // Get indices of relevant bin
                                        theta_index = get_index(theta,
                                                                max_theta,
                                                                num_theta);

                                        // Record particle's time on the surface
                                        A2_density[theta_index]
                                                  [phi__index_stuck]    +=
                                                time - t_stuck;
                                    }

                                    // Reset for this new bin
                                    t_stuck             = time;
                                    phi__index_stuck    = phi__index;
                                }
                            }
                        }
                    }

                }


                //========
                // Hop
                //========
                // Update theta_0, phi_0 for next hop
                theta_0 = theta;
                phi_0   = phi;

                // Reset position always to the eqm_source for at_altitude
                if (record_dist == "at_altitude") {
                    // Reset to keep phi_ unchanged (not phi)
                    theta_0 = theta_source;
                    phi_0   = phi_source - 2*pi*time/day;
                    if (phi_0 < 0) {
                        phi_0   += 2*pi;
                    }
                }

                // Generate random velocity magnitude and direction
                xi_1    = NR::ran3(seed);
                xi_2    = NR::ran3(seed);
                xi_3    = NR::ran3(seed);
                xi_4    = NR::ran3(seed);

                // Temperature and speed
                if (T_dist == "H_T_simple") {
                    // Find the phi_ time of day and temperature for that phi
                    phi_        = phi_prime(phi_0, time);
                    phi__index  = get_index(phi_, max_phi_, num_phi_);
                    theta_index = get_index(theta, max_theta, num_theta);

                    T           = A2_H_temperature[phi__index][theta_index];
                }
                else if (T_dist == "H_T_offsets") {
                    // First find and use the phi offset
                    phi_index   = get_index(phi_0, max_phi, num_phi);
                    theta_index = get_index(theta_0, max_theta, num_theta);

                    phi_offset  = A2_longitude_offset[phi_index][theta_index];

                    // Find the phi_ time and temperature for that offset phi
                    phi_        = phi_prime(phi_0 + phi_offset, time);
                    phi__index  = get_index(phi_, max_phi_, num_phi_);

                    T           = A2_H_temperature[phi__index][theta_index];
                }

                // Mare offset temperature
                if (misc == "mare_T_offset") {
                    if (is_in_mare(theta, phi, theta_m_1, theta_m_2,
                                   phi_m_1, phi_m_2)) {
                        T   += T_mare;
                    }

                    // Keep T below the maximum for generating v_0
                    if (T > 392) {
                        T   = 392;
                    }
                }

                v_0 = rand_vel(T, xi_1, xi_2, T_min, &A2_Maxwell_cdf);

                // Test for escape velocity, record and break if lost
                if (v_0 >= v_esc) {
                    N_escaped ++;
                    // Record when particle leaves
                    if (rates == "yes_rates") {
                        A1_time_esc.push_back(time);
                    }
                    else if (rates == "LADEE_lat_rates") {
                        if (is_LADEE_lat(theta)) {
                            A1_time_esc.push_back(time);
                        }
                    }

                    break;
                }

                // Generate direction angles
                eta = rand_eta(xi_3);

                if (emerge_dist == "M_B") {
                    alpha   = rand_alpha(xi_4);
                }
                else if (emerge_dist == "M_B_flux") {
                    alpha   = rand_alpha_flux(xi_4);
                }

                // Calculate flight variables
                a       = get_a(v_0);
                e       = get_e(a, alpha);
                h_max   = get_max_height(a, e);
                delta   = get_delta(a, e);
                t       = get_flight_time(a, e, delta);

                // Calculate end coordinates
                theta   = get_end_theta(theta_0, delta, eta);
                epsilon = get_epsilon(theta, theta_0, delta);
                phi     = get_end_phi(phi_0, epsilon, eta);

                // Add flight time
                time    += t;
                phi_    = phi_prime(phi, time);


                //========
                // Test for losses
                //========
                // Moving in/out of +/- 30 degrees latitude
                if (rates == "LADEE_lat_rates") {
                    // If was out and has now moved in
                    if ((not (is_LADEE_lat(theta_0))) &&
                        (is_LADEE_lat(theta))) {

                        A1_time_start.push_back(time);
                    }
                    // Or if was in and has now moved out
                    else if ((is_LADEE_lat(theta_0)) &&
                             (not (is_LADEE_lat(theta)))) {
                        A1_time_esc.push_back(time);
                    }
                }

                // Test for photodestruction (if on dayside) and capture
                if ((losses == "normal") || (losses == "normal&time")) {
                    // Photodestruction
                    if (is_in_day(phi_)) {
                        xi_1    = NR::ran3(seed);
                        Prob    = 1 - exp(-t/tau);
                        // Break if photodestroyed
                        if (xi_1 < Prob) {
                            N_destroyed ++;

                            // Record when particle leaves
                            if (rates == "yes_rates") {
                                A1_time_dest.push_back(time);
                            }
                            if (rates == "LADEE_lat_rates") {
                                if (is_LADEE_lat(theta)) {
                                    A1_time_dest.push_back(time);
                                }
                            }
                            break;
                        }
                    }

                    // PSR capture
                    if (cold_traps == "prob_at_theta") {
                        xi_1    = NR::ran3(seed);
                        if (test_capt_at_theta(theta, xi_1, f_N, f_S)) {
                            N_captured  ++;

                            // Record stuck on surface until end of simulation
                            if (record_dist == "density_stuck") {

                                // Store time stuck in the capture bin.
                                // Different notation: time is fixed so use \
                                    t_stuck as variable for changing bins.
                                // Particle leaves this first bin at time
                                t_stuck = day/(2*pi) * (1 - (phi_ -
                                                             floor(phi_)));

                                theta_index         = get_index(theta,
                                                                max_theta,
                                                                num_theta);
                                phi__index_stuck    = get_index(phi_, max_phi_,
                                                                num_phi_);

                                // Record time stuck
                                if (is_recording_time(time, N_bounce,
                                                      start_time,
                                                      start_N_bounce)) {
                                    A2_density[theta_index][phi__index_stuck]
                                            += t_stuck;
                                }

                                // Middle bins
                                // When simulation ends, will be at
                                phi__end    = phi + 2*pi*end_time/day;

                                // Sticks in each bin until the end from now or\
                                    start_time onwards, whichever is later
                                if (is_recording_time(time, N_bounce,
                                                      start_time,
                                                      start_N_bounce)) {
                                    for (int i_phi_=1; i_phi_<phi__end;
                                         i_phi_++) {
                                        phi__index  = phi__index_stuck + i_phi_;

                                        // Make between 0 and num_phi_
                                        phi__index  -= num_phi_ *
                                                       floor(phi__index /
                                                             num_phi_);

                                        // Record time stuck as Moon rotates
                                        A2_density[theta_index][phi__index] +=
                                                day/num_phi_;
                                    }
                                }
                                else {
                                    // Reset everything to be at start_time
                                    phi_    = phi_prime(phi, start_time);

                                    t_stuck = day/(2*pi) * (1 - (phi_ -
                                                                 floor(phi_)));

                                    phi__index_stuck    = get_index(phi_,
                                                                    max_phi_,
                                                                    num_phi_);

                                    // From start_time onwards
                                    for (int i_phi_=0; i_phi_<phi__end;
                                         i_phi_++) {
                                        phi__index  = phi__index_stuck + i_phi_;

                                        // Make between 0 and num_phi_
                                        phi__index  -= num_phi_ *
                                                       floor(phi__index /
                                                             num_phi_);

                                        // Record time stuck as Moon rotates
                                        A2_density[theta_index][phi__index] +=
                                                day/num_phi_;
                                    }
                                }


                                // Final bin
                                // Make between 0 and 2*pi
                                phi__end    -= 2*pi * floor(phi__end/(2*pi));

                                // Time spent in final bin:
                                t_stuck = day/(2*pi) * (phi__end -
                                                        floor(phi__end));

                                phi__index_stuck    = get_index(phi_, max_phi_,
                                                                num_phi_);

                                // Record time stuck
                                A2_density[theta_index][phi__index_stuck]   +=
                                        t_stuck;
                            }

                            // Record when particle leaves
                            if (rates == "yes_rates") {
                                A1_time_capt.push_back(time);
                            }

                            break;
                        }
                    }
                    else if (cold_traps == "local_traps") {
                        if (test_capt_local_traps(theta, phi, theta_trap_N_1,
                                                  theta_trap_N_2, phi_trap_N_1,
                                                  phi_trap_N_2, theta_trap_S_1,
                                                  theta_trap_S_2, phi_trap_S_1,
                                                  phi_trap_S_2)) {
                            N_captured  ++;

                            // Record when particle leaves
                            if (rates == "yes_rates") {
                                A1_time_capt.push_back(time);
                            }

                            break;
                        }
                    }
                    else if (cold_traps == "seasonal") {
                        // If captured in a permanent trap, then normal loss
                        xi_1    = NR::ran3(seed);
                        if (test_capt_at_theta(theta, xi_1, f_N, f_S)) {
                            N_captured  ++;

                            // Record when particle leaves
                            if (rates == "yes_rates") {
                                A1_time_capt.push_back(time);
                            }

                            break;
                        }

                        // If captured in a seasonal trap, then reappear later
                        xi_2    = NR::ran3(seed);
                        if (test_capt_seasonal(theta, xi_2, time, f_season_N,
                                               f_season_S)) {
                            // Record time of capture
                            if ((rates == "yes_rates") ||
                                (rates == "LADEE_lat_rates")) {

                                if (is_north_pole(theta)) {
                                    A1_time_seasonal_capt_N.push_back(time);
                                }
                                else if (is_south_pole(theta)) {
                                    A1_time_seasonal_capt_S.push_back(time);
                                }
                            }

                            // Release when the trap shrinks
                            xi_3    = NR::ran3(seed);
                            time    = get_time_seasonal_rel(theta, time, xi_3);

                            // Record time of release
                            if ((rates == "yes_rates") ||
                                (rates == "LADEE_lat_rates")) {

                                if (is_north_pole(theta)) {
                                    A1_time_seasonal_rel_N.push_back(time);
                                }
                                else if (is_south_pole(theta)) {
                                    A1_time_seasonal_rel_S.push_back(time);
                                }
                            }
                        }
                    }
                }

                // Solar flare loss event (destroys all dayside particles)
                if (misc == "solar_flare") {
                    // If not at night and during the time of the flare
                    if ((is_in_day(phi_)) && (t_flare_1 < time) &&
                        (time < t_flare_2)) {

                        N_destroyed ++;

                        // Record when particle leaves
                        if (rates == "yes_rates") {
                            A1_time_dest.push_back(time);
                        }
                        break;
                    }
                }


                //========
                // Record particle density
                //========
                // Add the particle's contribution to the density distribution
                if (record_dist == "density") {
                    if (is_recording_time(time, N_bounce, start_time,
                                          start_N_bounce)) {
                        // Get indices of relevant bin in the array
                        theta_index = get_index(theta, max_theta, num_theta);
                        phi__index  = get_index(phi_, max_phi_, num_phi_);

                        // Record flux of the particle as it reaches the surface
                        flux    = get_flux(v_0, alpha);
                        A2_density[theta_index][phi__index] += flux;

                        // Record also from when the particle left the surface
                        theta_index = get_index(theta_0, max_theta, num_theta);
                        phi__1      = phi_prime(phi_0, time-t);
                        phi__index  = get_index(phi__1, max_phi_, num_phi_);

                        A2_density[theta_index][phi__index] += flux;
                        // For making an eqm_source distribution, use only
//                        A2_density[theta_index][phi__index] += 1;
                    }
                }

                else if (record_dist == "densities") {
                    if (is_recording_time(time, N_bounce, start_time,
                                          start_N_bounce)) {
                        // Get time of day index to select which density array
                        phi__index  = get_index(phi_, max_phi_, num_phi_);
                        // Get indices of relevant bin in the array
                        theta_index = get_index(theta, max_theta, num_theta);
                        phi_index   = get_index(phi, max_phi, num_phi);

                        // Record flux of the particle as it reaches the surface
                        flux    = get_flux(v_0, alpha);
                        A3_density[phi__index][theta_index][phi_index]  += flux;

                        // Record also from when the particle left the surface
                        phi__1      = phi_prime(phi_0, time - t);
                        phi__index  = get_index(phi__1, max_phi_, num_phi_);

                        theta_index = get_index(theta_0, max_theta, num_theta);
                        phi_index   = get_index(phi_0, max_phi, num_phi);

                        A3_density[phi__index][theta_index][phi_index]  += flux;
                    }
                }

                else if (record_dist == "all_densities") {
                    if (is_recording_time(time, N_bounce, start_time,
                                          start_N_bounce)) {
                        // Get time index to select which density array
                        time_index  = get_index(time - start_time,
                                                num_day * day,
                                                num_day * num_time);
                        // Get indices of relevant bin in the array
                        theta_index = get_index(theta, max_theta, num_theta);
                        phi__index  = get_index(phi_, max_phi_, num_phi_);

                        // Record flux of the particle as it reaches the surface
                        flux    = get_flux(v_0, alpha);
                        A3_density[time_index][theta_index][phi__index] += flux;

                        // Record also from when the particle left the surface
                        time_index  = get_index(time - t - start_time,
                                                num_day * day,
                                                num_day * num_time);
                        if (time_index < 0) {
                            time_index  = 0;
                        }

                        theta_index = get_index(theta_0, max_theta, num_theta);
                        phi__1      = phi_prime(phi_0, time-t);
                        phi__index  = get_index(phi__1, max_phi_, num_phi_);

                        A3_density[time_index][theta_index][phi__index] += flux;
                    }
                }

                else if (record_dist == "at_altitude") {
                    if (is_recording_time(time, N_bounce, start_time,
                                          start_N_bounce)) {

                        // Find altitude step of the peak
                        alt_step_max    = floor(h_max / alt_step);

                        // Density at the surface (i_alt = 0)
                        // Before peak (launching)
                        theta_index = get_index(theta_0, max_theta, num_theta);
                        phi__1      = phi_prime(phi_0, time-t);
                        phi__index  = get_index(phi__1, max_phi_, num_phi_);

                        // Record flux of the particle as it reaches the surface
                        flux    = get_flux(v_0, alpha);

                        A4_density_alt[theta_index][phi__index][0][i_run]   +=
                                flux;

                        // After peak (landing)
                        theta_index = get_index(theta, max_theta, num_theta);
                        phi__index  = get_index(phi_, max_phi_, num_phi_);

                        A4_density_alt[theta_index][phi__index][0][i_run]   +=
                                flux;


                        // First altitude step (i_alt = 1)
                        if (alt_step_max > 0) {
                            alt     = alt_step;

                            // Time of flight wrt peak
                            delta_p = get_delta_p(a, e, alt);
                            t_p     = get_flight_time(a, e, delta_p);

                            // Calculate the speed at the step altitude
                            v       = get_v(a, alt);
                            // and the angle away from horizontal
                            alpha   = get_alpha(a, e, alt);

                            // Before peak
                            delta_1 = delta/2 - delta_p;
                            theta_1 = get_end_theta(theta_0, delta_1, eta);
                            epsilon = get_epsilon(theta_1, theta_0,  delta_1);
                            phi_1   = get_end_phi(phi_0, epsilon, eta);
                            phi__1  = phi_prime(phi_1, time-(0.5*t+t_p));

                            theta_index = get_index(theta_1, max_theta,
                                                    num_theta);
                            phi__index  = get_index(phi__1, max_phi_, num_phi_);

                            // Weight by time spent near this altitude
                            flux    = get_flux(v, alpha);

                            // Store values
                            A4_density_alt[theta_index][phi__index][1][i_run]
                                    += flux;

                            // After peak
                            delta_2 = delta - delta_1;
                            theta_2 = get_end_theta(theta_0, delta_2, eta);
                            epsilon = get_epsilon(theta_2, theta_0, delta_2);
                            phi_2   = get_end_phi(phi_0, epsilon, eta);
                            phi__2  = phi_prime(phi_2, time-(0.5*t-t_p));

                            theta_index = get_index(theta_2, max_theta,
                                                    num_theta);
                            phi__index  = get_index(phi__2, max_phi_, num_phi_);

                            A4_density_alt[theta_index][phi__index][1][i_run]
                                    += flux;
                        }

                        // All other altitude steps that the particle crosses
                        int i_alt = 2;
                        while ((i_alt <= alt_step_max) && (i_alt < num_alt)) {
                            alt     = i_alt * alt_step;

                            // delta_p at the step altitude
                            delta_p = get_delta_p(a, e, alt);
                            t_p     = get_flight_time(a, e, delta_p);

                            // Calculate the speed at the step altitude
                            v       = get_v(a, alt);
                            // and the angle away from horizontal
                            alpha   = get_alpha(a, e, alt);

                            // Before peak
                            delta_1 = delta/2 - delta_p;
                            theta_1 = get_end_theta(theta_0, delta_1, eta);
                            epsilon = get_epsilon(theta_1, theta_0,  delta_1);
                            phi_1   = get_end_phi(phi_0, epsilon, eta);
                            phi__1  = phi_prime(phi_1, time-(0.5*t+t_p));

                            theta_index = get_index(theta_1, max_theta,
                                                    num_theta);
                            phi__index  = get_index(phi__1, max_phi_, num_phi_);

                            // Weight by time spent near this altitude
                            flux    = get_flux(v, alpha);

                            // Store
                            A4_density_alt[theta_index][phi__index][i_alt]
                                          [i_run]   += flux;

                            // After peak
                            delta_2 = delta - delta_1;
                            theta_2 = get_end_theta(theta_0, delta_2, eta);
                            epsilon = get_epsilon(theta_2, theta_0, delta_2);
                            phi_2   = get_end_phi(phi_0, epsilon, eta);
                            phi__2  = phi_prime(phi_2, time-(0.5*t-t_p));

                            theta_index = get_index(theta_2, max_theta,
                                                    num_theta);
                            phi__index  = get_index(phi__2, max_phi_, num_phi_);

                            A4_density_alt[theta_index][phi__index][i_alt]
                                          [i_run]   += flux;

                            // Next altitude step
                            i_alt   ++;
                        }
                    }
                }


                // Test for end of simulation by time or number of bounces
                if ((losses == "time") || (losses == "normal&time")) {
                    if ((N_bounce > end_N_bounce) && (time > end_time)) {
                        N_ended ++;
                        break;
                    }
                }
            }
            //                              ================================= //
            //                                              Bounce until loss //
            //                              ================================= //
        }
        //                                  ================================= //
        //                                                      Particle Runs //
        //                                  ================================= //

        // Store fraction of captured/destroyed/escaped particles from each run
        frac_captured   = (double) N_captured/num_particle;
        frac_destroyed  = (double) N_destroyed/num_particle;
        frac_escaped    = (double) N_escaped/num_particle;
        A1_frac_captured[i_run]     = frac_captured;
        A1_frac_destroyed[i_run]    = frac_destroyed;
        A1_frac_escaped[i_run]      = frac_escaped;
    }
    //                                      ================================= //
    //                                                        Simulation Runs //
    //                                      ================================= //


    // Print when finished
    cout << "[" << world_rank << "]: " << runs_per_proc << " runs completed."
         << endl;


    //========
    // Results
    //========
    // Calculate loss fractions
    mean_frac_captured  = get_mean(A1_frac_captured);
    mean_frac_destroyed = get_mean(A1_frac_destroyed);
    mean_frac_escaped   = get_mean(A1_frac_escaped);
    stdev_captured      = get_stdev(A1_frac_captured);
    stdev_destroyed     = get_stdev(A1_frac_destroyed);
    stdev_escaped       = get_stdev(A1_frac_escaped);

    // Make list of all means and vars=stdevs^2 for sending and receiving
    double  A1_mean_var[6]  = {mean_frac_captured, mean_frac_destroyed,
                               mean_frac_escaped, stdev_captured*stdev_captured,
                               stdev_destroyed*stdev_destroyed,
                               stdev_escaped*stdev_escaped};


    // Altitude density means and stdevs
    if (record_dist == "at_altitude") {
        double stdev_tmp;
        for (int i_theta=0; i_theta<num_theta; i_theta++) {
            for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                for (int i_alt=0; i_alt<num_alt; i_alt++) {
                    A3_density_alt[i_theta][i_phi_][i_alt]    = get_mean(
                            A4_density_alt[i_theta][i_phi_][i_alt]);

                    stdev_tmp   = get_stdev(A4_density_alt[i_theta][i_phi_]
                                                          [i_alt]);

                    A3_variance_alt[i_theta][i_phi_][i_alt] = stdev_tmp *
                                                              stdev_tmp;
                }
            }
        }
    }


    //=========================
    // [1], ..., [N]: Send results to [0]
    //=========================
    if (world_rank > 0) {
        // Send means and stdevs
        MPI_Send(A1_mean_var, 6, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);

        // Send density vectors (one at a time with own tag)
        if ((record_dist == "density") || (record_dist == "density_stuck")) {
            for (int i_theta=0; i_theta<num_theta; i_theta++) {
                MPI_Send(&A2_density[i_theta].front(), num_phi_, MPI_DOUBLE, 0,
                         10 + i_theta, MPI_COMM_WORLD);
            }
        }

        // Send eqm_density vectors
        else if (record_dist == "densities") {
            // For each density array
            for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                // Send density vectors (one at a time with own tag)
                for (int i_theta=0; i_theta<num_theta; i_theta++) {
                    MPI_Ssend(&A3_density[i_phi_][i_theta].front(), num_phi,
                              MPI_DOUBLE, 0, 10 + i_theta, MPI_COMM_WORLD);
                }
            }
        }

        // Send time order density vectors
        else if (record_dist == "all_densities") {
            // For each density array
            for (int i_day=0; i_day<num_day; i_day++) {
                for (int i_t=0; i_t<num_time; i_t++) {
                    // Send density vectors (one at a time with own tag)
                    for (int i_theta=0; i_theta<num_theta; i_theta++) {
                        MPI_Ssend(&A3_density[(i_day*num_time) + i_t]
                                             [i_theta].front(),
                                  num_phi_, MPI_DOUBLE, 0, 10 + i_theta,
                                  MPI_COMM_WORLD);
                    }

                }
            }
        }

        // Send altitude density vectors
        else if (record_dist == "at_altitude") {
            // For each density array
            for (int i_theta=0; i_theta<num_theta; i_theta++) {
                for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                    // Send alt density vectors (one at a time with own tag)
                    MPI_Ssend(&A3_density_alt[i_theta][i_phi_].front(),
                              num_alt, MPI_DOUBLE, 0, 10 + i_phi_,
                              MPI_COMM_WORLD);
                }
            }
            // For each density variance array
            for (int i_theta=0; i_theta<num_theta; i_theta++) {
                for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                    // Send alt variance vectors (one at a time with own tag)
                    MPI_Ssend(&A3_variance_alt[i_theta][i_phi_].front(),
                              num_alt, MPI_DOUBLE, 0, 10 + i_phi_,
                              MPI_COMM_WORLD);
                }
            }
        }


        // Send start and end times etc.
        if ((rates == "yes_rates") || (rates == "LADEE_lat_rates")) {
            MPI_Send(&A1_time_start.front(), A1_time_start.size(),
                     MPI_DOUBLE, 0, 5, MPI_COMM_WORLD);
            MPI_Send(&A1_time_capt.front(), A1_time_capt.size(),
                     MPI_DOUBLE, 0, 6, MPI_COMM_WORLD);
            MPI_Send(&A1_time_dest.front(), A1_time_dest.size(),
                     MPI_DOUBLE, 0, 7, MPI_COMM_WORLD);
            MPI_Send(&A1_time_esc.front(), A1_time_esc.size(),
                     MPI_DOUBLE, 0, 8, MPI_COMM_WORLD);
            MPI_Send(&A1_time_seasonal_capt_N.front(),
                     A1_time_seasonal_capt_N.size(),
                     MPI_DOUBLE, 0, 9, MPI_COMM_WORLD);
            MPI_Send(&A1_time_seasonal_rel_N.front(),
                     A1_time_seasonal_rel_N.size(),
                     MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
            MPI_Send(&A1_time_seasonal_capt_S.front(),
                     A1_time_seasonal_capt_S.size(),
                     MPI_DOUBLE, 0, 11, MPI_COMM_WORLD);
            MPI_Send(&A1_time_seasonal_rel_S.front(),
                     A1_time_seasonal_rel_S.size(),
                     MPI_DOUBLE, 0, 12, MPI_COMM_WORLD);
        }
    }


    //=========================
    // [0]: Collect and process results from [1], ..., [N]
    //      Save results
    //=========================
    else if (world_rank == 0) {
        int             buffer_length;
        vector<double>  A1_buffer;

        MPI_Status  status;


        //========
        // Collect from each other processor
        //========
        // Inefficient to loop over all ranks in order but it's simple and \
           is negligible compared to the processing time.
        for (int i_rank=1; i_rank<world_size; i_rank++) {
            // Receive means and vars
            A1_buffer.resize(6);
            MPI_Recv(&A1_buffer.front(), 6, MPI_DOUBLE, i_rank, 0,
                     MPI_COMM_WORLD, MPI_STATUS_IGNORE);

            // Add to the totals
            for (int i_var=0; i_var<6; i_var++) {
                A1_mean_var[i_var]  += A1_buffer[i_var];
            }

            if ((record_dist == "density") ||
                (record_dist == "density_stuck")) {

                A1_buffer.resize(num_phi_);
                // Receive density vectors and add to total density
                for (int i_theta=0; i_theta<num_theta; i_theta++) {
                    MPI_Recv(&A1_buffer.front(), num_phi_, MPI_DOUBLE,
                             i_rank, 10 + i_theta, MPI_COMM_WORLD,
                             MPI_STATUS_IGNORE);

                    A2_density[i_theta] += A1_buffer;
                }
            }
            else if (record_dist == "densities") {
                A1_buffer.resize(num_phi);
                // For each density array
                for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                    // Receive density vectors and add to total density
                    for (int i_theta=0; i_theta<num_theta; i_theta++) {
                        MPI_Recv(&A1_buffer.front(), num_phi, MPI_DOUBLE,
                                 i_rank, 10 + i_theta, MPI_COMM_WORLD,
                                 MPI_STATUS_IGNORE);

                        A3_density[i_phi_][i_theta] += A1_buffer;
                    }
                }
            }
            else if (record_dist == "all_densities") {
                A1_buffer.resize(num_phi_);
                // For each density array
                for (int i_day=0; i_day<num_day; i_day++) {
                    for (int i_t=0; i_t<num_time; i_t++) {
                        // Receive density vectors and add to total density
                        for (int i_theta=0; i_theta<num_theta; i_theta++) {
                            MPI_Recv(&A1_buffer.front(), num_phi_,
                                     MPI_DOUBLE, i_rank, 10 + i_theta,
                                     MPI_COMM_WORLD, MPI_STATUS_IGNORE);

                            A3_density[(i_day*num_time) + i_t][i_theta] +=
                                    A1_buffer;
                        }
                    }
                }
            }
            else if (record_dist == "at_altitude") {
                A1_buffer.resize(num_alt);
                // For each density bin
                for (int i_theta=0; i_theta<num_theta; i_theta++) {
                    for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                        // Receive vectors and add to total density
                        MPI_Recv(&A1_buffer.front(), num_alt, MPI_DOUBLE,
                                 i_rank, 10 + i_phi_, MPI_COMM_WORLD,
                                 MPI_STATUS_IGNORE);

                        A3_density_alt[i_theta][i_phi_] += A1_buffer;
                    }
                }
                // For each density variances bin
                for (int i_theta=0; i_theta<num_theta; i_theta++) {
                    for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                        // Receive vectors and add to total density
                        MPI_Recv(&A1_buffer.front(), num_alt, MPI_DOUBLE,
                                 i_rank, 10 + i_phi_, MPI_COMM_WORLD,
                                 MPI_STATUS_IGNORE);

                        A3_variance_alt[i_theta][i_phi_]    += A1_buffer;
                    }
                }
            }

            // Receive start times etc.
            if ((rates == "yes_rates") || (rates == "LADEE_lat_rates")) {
                // Unknown lengths!
                // Probe for a sent vector and set the required buffer length
                MPI_Probe(i_rank, 5, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                // Receive vector and append it to rank 0's vector
                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         5, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                // Append to rank 0's vector
                A1_time_start.insert(A1_time_start.end(), A1_buffer.begin(),
                                     A1_buffer.end());

                // Repeat for capture
                MPI_Probe(i_rank, 6, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         6, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                A1_time_capt.insert(A1_time_capt.end(), A1_buffer.begin(),
                                    A1_buffer.end());

                // Repeat for destruction
                MPI_Probe(i_rank, 7, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         7, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                A1_time_dest.insert(A1_time_dest.end(), A1_buffer.begin(),
                                    A1_buffer.end());

                // Repeat for escape
                MPI_Probe(i_rank, 8, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         8, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                A1_time_esc.insert(A1_time_esc.end(), A1_buffer.begin(),
                                   A1_buffer.end());

                // Repeat for seasonal N capture
                MPI_Probe(i_rank, 9, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         9, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                A1_time_seasonal_capt_N.insert(A1_time_seasonal_capt_N.end(),
                                               A1_buffer.begin(),
                                               A1_buffer.end());

                // Repeat for seasonal N release
                MPI_Probe(i_rank, 10, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                A1_time_seasonal_rel_N.insert(A1_time_seasonal_rel_N.end(),
                                              A1_buffer.begin(),
                                              A1_buffer.end());

                // Repeat for seasonal S capture
                MPI_Probe(i_rank, 11, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         11, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                A1_time_seasonal_capt_S.insert(A1_time_seasonal_capt_S.end(),
                                               A1_buffer.begin(),
                                               A1_buffer.end());

                // Repeat for seasonal S release
                MPI_Probe(i_rank, 12, MPI_COMM_WORLD, &status);
                MPI_Get_count(&status, MPI_DOUBLE, &buffer_length);
                A1_buffer.resize(buffer_length);

                MPI_Recv(&A1_buffer.front(), buffer_length, MPI_DOUBLE, i_rank,
                         12, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                A1_time_seasonal_rel_S.insert(A1_time_seasonal_rel_S.end(),
                                              A1_buffer.begin(),
                                              A1_buffer.end());
            }

            cout << "[0]: Recvd from [" << i_rank << "]" << endl;
        }


        // Divide densities by recording time to convert to num per bin per m, \
            and by the source rate so can later be converted into real particles
        if ((record_dist == "density") || (record_dist == "density_stuck")) {
            for (int i_theta=0; i_theta<num_theta; i_theta++) {
                for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                    A2_density[i_theta][i_phi_] /= end_time - start_time;
                    A2_density[i_theta][i_phi_] /= N_sim;
                }
            }
        }
        else if (record_dist == "densities") {
            for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                for (int i_theta=0; i_theta<num_theta; i_theta++) {
                    for (int i_phi=0; i_phi<num_phi; i_phi++) {
                        A3_density[i_phi_][i_theta][i_phi]  /= end_time -
                                                               start_time;
                        A3_density[i_phi_][i_theta][i_phi]  /= N_sim;
                    }
                }
            }
        }
        else if (record_dist == "all_densities") {
            for (int i_day=0; i_day<num_day; i_day++) {
                for (int i_time=0; i_time<num_time; i_time++) {
                    for (int i_theta=0; i_theta<num_theta; i_theta++) {
                        for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                            A3_density[(i_day*num_time) + i_time][i_theta]
                                      [i_phi_]  /= day / num_time;
                            A3_density[(i_day*num_time) + i_time][i_theta]
                                      [i_phi_]  /= N_sim;
                        }
                    }
                }
            }
        }
        else if (record_dist == "at_altitude") {
            // Densities and variances
            for (int i_theta=0; i_theta<num_theta; i_theta++) {
                for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                    for (int i_alt=0; i_alt<num_alt; i_alt++) {
                        A3_density_alt[i_theta][i_phi_][i_alt]  /=
                                end_time - start_time;
                        A3_variance_alt[i_theta][i_phi_][i_alt] /=
                                (end_time - start_time)*(end_time - start_time);

                        A3_density_alt[i_theta][i_phi_][i_alt]  /= N_sim;
                        A3_variance_alt[i_theta][i_phi_][i_alt] /= N_sim*N_sim;
                    }
                }
            }
        }


        // Calculate final loss results
        for (int i_var=0; i_var<16; i_var++) {
            A1_mean_var[i_var]  /= world_size;
        }

        // Unpack means and stdevs (back from vars) from array
        mean_frac_captured  = A1_mean_var[0];
        mean_frac_destroyed = A1_mean_var[1];
        mean_frac_escaped   = A1_mean_var[2];
        stdev_captured      = sqrt(A1_mean_var[3]);
        stdev_destroyed     = sqrt(A1_mean_var[4]);
        stdev_escaped       = sqrt(A1_mean_var[5]);

        // Sort the "yes_rates"/"LADEE_lat_rates" times into bins to get rates
        if (end_time > 0) {
            A1_bin_time_start   = bin_values(A1_time_start, num_time_rates,
                                             end_time);
            A1_bin_time_capt    = bin_values(A1_time_capt, num_time_rates,
                                             end_time);
            A1_bin_time_dest    = bin_values(A1_time_dest, num_time_rates,
                                             end_time);
            A1_bin_time_esc     = bin_values(A1_time_esc, num_time_rates,
                                             end_time);
            A1_bin_time_seasonal_capt_N = bin_values(A1_time_seasonal_capt_N,
                                                     num_time_rates, end_time);
            A1_bin_time_seasonal_rel_N  = bin_values(A1_time_seasonal_rel_N,
                                                     num_time_rates, end_time);
            A1_bin_time_seasonal_capt_S = bin_values(A1_time_seasonal_capt_S,
                                                     num_time_rates, end_time);
            A1_bin_time_seasonal_rel_S  = bin_values(A1_time_seasonal_rel_S,
                                                     num_time_rates, end_time);

            // Print time size of each bin for units
            if (world_rank == 0) {
                cout << endl << "Rates bin size = " << end_time / num_time_rates
                     << " s" << endl;
            }
        }


        //========
        // Print and save results
        //========
        // Loss fractions
        cout << endl;
        print_debug6(mean_frac_captured, stdev_captured,
                     mean_frac_destroyed, stdev_destroyed,
                     mean_frac_escaped, stdev_escaped);
        cout << endl;

        cout << "Saving files..." << endl << endl;

        // Save results directly to the (current) results folder
        string  dir_cur = dir;
        dir_cur         = dir_cur.append("results/");

        // Density distributions
        if ((record_dist == "density") || (record_dist == "density_stuck")) {
            // Save normally if not part of an array of jobs
            if (job_index == 0) {
                save_to_dir_vv("density", dir_cur, 'n', A2_density);
            }
            // Otherwise, save with the job index label
            else {
                string  file_with_index = "density_";

                file_with_index = file_with_index.append(SSTR(job_index));

                save_to_dir_vv(file_with_index, dir_cur, 'n', A2_density);
            }
        }
        else if (record_dist == "densities") {
            string  dir_den = dir_cur;
            dir_den         = dir_den.append("densities/");
            string  filename;
            // Save each density map to its own file
            for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                filename    = "density_";
                filename.append(SSTR(i_phi_));
                save_to_dir_vv(filename.c_str(), dir_den, 'n',
                               A3_density[i_phi_]);
            }
        }
        else if (record_dist == "all_densities") {
            string  dir_den = dir_cur;
            dir_den         = dir_den.append("all_densities/");
            string  filename;
            // Save each density map to its own file
            for (int i_day=0; i_day<num_day; i_day++) {
                for (int i_t=0; i_t<num_time; i_t++) {
                    filename    = "density_";
                    filename.append(SSTR(i_day));
                    filename.append("_");
                    filename.append(SSTR(i_t));
                    save_to_dir_vv(filename.c_str(), dir_den, 'n',
                                   A3_density[(i_day*num_time) + i_t]);
                }
            }
        }
        else if (record_dist == "at_altitude") {
            // No job_index
//            save_to_dir_vvv("alt_density", dir_cur, 'n', A3_density_alt);
//            save_to_dir_vvv("variance_alt", dir_cur, 'n', A3_variance_alt);

            // Multi run saving using the job_index
            string  file_with_index     = "alt_density/alt_density_";
            string  file_with_index2    = "alt_density/alt_vars_";

            file_with_index     = file_with_index.append(SSTR(job_index));
            file_with_index2    = file_with_index2.append(SSTR(job_index));

            save_to_dir_vvv(file_with_index, dir_cur, 'n', A3_density_alt);
            save_to_dir_vvv(file_with_index2, dir_cur, 'n', A3_variance_alt);
        }

        // Start and end times
        if ((rates == "yes_rates") || (rates == "LADEE_lat_rates")) {
            save_to_dir_v("bin_time_start", dir_cur, 'n', A1_bin_time_start);
            save_to_dir_v("bin_time_capt", dir_cur, 'n', A1_bin_time_capt);
            save_to_dir_v("bin_time_dest", dir_cur, 'n', A1_bin_time_dest);
            save_to_dir_v("bin_time_esc", dir_cur, 'n', A1_bin_time_esc);
            save_to_dir_v("bin_time_seasonal_capt_N", dir_cur, 'n',
                          A1_bin_time_seasonal_capt_N);
            save_to_dir_v("bin_time_seasonal_rel_N", dir_cur, 'n',
                          A1_bin_time_seasonal_rel_N);
            save_to_dir_v("bin_time_seasonal_capt_S", dir_cur, 'n',
                          A1_bin_time_seasonal_capt_S);
            save_to_dir_v("bin_time_seasonal_rel_S", dir_cur, 'n',
                          A1_bin_time_seasonal_rel_S);
        }


        //========
        // Copy or save everything to a new folder in the results_backup folder.
        //========
        if (backup_results == "yes_backup") {

            cout << "Saving backup files..." << endl << endl;
            // Make the new folder named by job ID and current date and time
            string  dir_bak = dir;
            dir_bak.append("results_backup/");
            dir_bak.append(SSTR(job_ID));
            if (job_index != 0) {
                dir_bak.append("_");
                dir_bak.append(SSTR(job_index));
            }
            mkdir(dir_bak.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

            dir_bak = dir_bak.append("/");

            // Make new sub-folder for the densities
            string  dir_bak_den = dir_bak;
            string  dir_bak_alt = dir_bak;
            dir_bak_den.append("densities/");
            dir_bak_alt.append("alt_density/");
            mkdir(dir_bak_den.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
            mkdir(dir_bak_alt.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

            // Save distributions directly to the new folder
            if ((record_dist == "density") ||
                (record_dist == "density_stuck")) {
                save_to_dir_vv("density", dir_bak, 'n', A2_density);
            }
            else if (record_dist == "densities") {
                string  filename;
                // Save each density map to its own file
                for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
                    filename    = "density_";
                    filename.append(SSTR(i_phi_));
                    save_to_dir_vv(filename.c_str(), dir_bak_den, 'n',
                                   A3_density[i_phi_]);
                }
            }
            else if (record_dist == "all_densities") {
                string  filename;
                // Save each density map to its own file
                for (int i_day=0; i_day<num_day; i_day++) {
                    for (int i_t=0; i_t<num_time; i_t++) {
                        filename    = "density_";
                        filename.append(SSTR(i_day));
                        filename.append("_");
                        filename.append(SSTR(i_t));
                        save_to_dir_vv(filename.c_str(), dir_bak_den, 'n',
                                       A3_density[(i_day*num_time) + i_t]);
                    }
                }
            }
            else if (record_dist == "at_altitude") {
                // No job_index
//                save_to_dir_vvv("alt_density", dir_bak, 'n', A3_density_alt);
//                save_to_dir_vvv("variance_alt", dir_bak, 'n', A3_variance_alt);

                // Multi run saving using the job_index
                string  file_with_index     = "alt_density/alt_density_";
                string  file_with_index2    = "alt_density/alt_vars_";

                file_with_index     = file_with_index.append(SSTR(job_index));
                file_with_index2    = file_with_index2.append(SSTR(job_index));

                save_to_dir_vvv(file_with_index, dir_bak, 'n', A3_density_alt);
                save_to_dir_vvv(file_with_index2, dir_bak, 'n',
                                A3_variance_alt);
            }

            // Save start and end times directly to the new folder
            if ((rates == "yes_rates") || (rates == "LADEE_lat_rates")) {
                save_to_dir_v("bin_time_start", dir_bak, 'n',
                              A1_bin_time_start);
                save_to_dir_v("bin_time_capt", dir_bak, 'n', A1_bin_time_capt);
                save_to_dir_v("bin_time_dest", dir_bak, 'n', A1_bin_time_dest);
                save_to_dir_v("bin_time_esc", dir_bak, 'n', A1_bin_time_esc);
                save_to_dir_v("bin_time_seasonal_capt_N", dir_bak, 'n',
                              A1_bin_time_seasonal_capt_N);
                save_to_dir_v("bin_time_seasonal_rel_N", dir_bak, 'n',
                              A1_bin_time_seasonal_rel_N);
                save_to_dir_v("bin_time_seasonal_capt_S", dir_bak, 'n',
                              A1_bin_time_seasonal_capt_S);
                save_to_dir_v("bin_time_seasonal_rel_S", dir_bak, 'n',
                              A1_bin_time_seasonal_rel_S);
            }
        }


        // Calculate the total time elapsed
        t_end   = clock();
        double  t_elapsed   = (double) (t_end - t_start) / CLOCKS_PER_SEC;
        cout << endl;
        print_debug(t_elapsed);
        cout << "Elapsed time " << floor(t_elapsed/(60*60)) << ':'
             << floor(t_elapsed/60 - floor(t_elapsed/(60*60)) * 60) << ':'
             << floor(t_elapsed - floor(t_elapsed/(60)) * 60) << endl;

        // Print the end time
        cout << "Ended at ";
        print_date_time('y', 'y');
        cout << endl << endl;
    }


    // ======================== End MPI
    MPI_Finalize();
    // ================================================================== //
}
//=========






























