/**
    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

    Functions.
**/


/*
    0. Includes etc.
    1. Constants
    2. Functions
        2.1 RNG
        2.2 Migration
        2.3 Cold Traps
        2.4 Misc
        2.5 Admin etc.
*/


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


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

#include <iostream>
#include <cmath>
#include <numeric>
#include <ctime>

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. Functions                                  //
////////////////////////////////////////////////////////////////////////////////


//============================================================================//
//                              2.1 RNG                                       //
//============================================================================//


//=========================
// Sources
//=========================
double  rand_phi_iso(double xi)
{
    /// Returns a random isotropic angle phi from a random number 0 <= xi < 1 \
        using the inverse cdf for phi.
    return  2*pi*xi;
}
//========


double  rand_theta_iso(double xi)
{
    /// Returns a random isotropic angle theta from a random number 0 <= xi < 1\
        using the inverse cdf for theta.
    return  acos(1 - 2*xi);
}
//========


double rand_phi_solar(double xi_1, double xi_2, double xi_3)
{
    /// Returns a random angle phi from two random numbers \
        using the Crider and Vodrak 2000 solar wind source distribution.
    double  phi;

    phi = asin(sqrt((xi_1 * cos(2*pi * xi_2) * cos(2*pi * xi_2)) /
                    (1 - xi_1 * sin(2*pi * xi_2) * sin(2*pi * xi_2))));

    // Even chance either side of subsolar point
    if (xi_3 < 0.5) {
        phi = 2*pi - phi;
    }

    return  phi;
}
//========


double  rand_theta_solar(double xi_1, double xi_2)
{
    /// Returns a random angle theta from two random numbers \
        using the Crider and Vodrak 2000 solar wind source distribution.
    return  acos(sqrt(xi_1) * sin(2*pi * xi_2));
}
//========


//=========================
// Velocity
//=========================
double  Maxwell_cdf(int T, double v, const double m)
{
    /// Calculates the Maxwell cdf for a mass m, Temperature T, and velocity v.
    double  gamma   = m/(2*k*T);

    return  erf(sqrt(gamma) * v) - (2 / sqrt(pi)) * sqrt(gamma) * v *
            exp(-gamma * v*v);
}
//========


double  Maxwell_flux_cdf(int T, double v, const double m)
{
    /// Calculates the Maxwell flux cdf for a mass m, Temperature T, velocity v.
    double  gamma   = m/(2*k*T);

    return  1 - (1 + gamma * v*v) * exp(-gamma * v*v);
}
//========


vector< vector<double> >    make_A2_Maxwell_cdf(double T_0, double T_1,
                                                 double T_min, const double m)
{
    /// Returns a 2D array of Maxwell cdf values for different T, v: \
        T from T_min to T_min+T_0+T_1 in 1K increments. \
        v from 0 to v_esc in v_esc/1000 increments.
    vector< vector<double> >    A2_Maxwell_cdf;
    int     T,
            // Number of temperatures (inclusive)
            num_T   = T_0 + T_1 - T_min + 1;
    double  v;

    A2_Maxwell_cdf.resize(num_T);
    // Loop over each temperature
    for (int i_T=0; i_T<num_T; i_T++) {
        // T goes up in 1K increments
        T   = i_T + T_min;
        // Loop over each velocity
        for (int i_v=0; i_v<=1001; i_v++) {
            // v up in v_esc/1000 increments
            v   = i_v * v_esc/1000;
            A2_Maxwell_cdf[i_T].push_back(Maxwell_cdf(T, v, m));
        }
    }

    return  A2_Maxwell_cdf;
}
//========


vector< vector<double> >    make_A2_Maxwell_flux_cdf(double T_0, double T_1,
                                                      double T_min,
                                                      const double m)
{
    /// Returns a 2D array of Maxwell flux cdf values for different T, v: \
        T from T_min to T_min+T_0+T_1 in 1K increments. \
        v from 0 to v_esc in v_esc/1000 increments.
    vector< vector<double> >    A2_Maxwell_cdf;
    int     T,
            // Number of temperatures (inclusive)
            num_T   = T_0 + T_1 - T_min + 1;
    double  v;

    A2_Maxwell_cdf.resize(num_T);
    // Loop over each temperature
    for (int i_T=0; i_T<num_T; i_T++) {
        // T goes up in 1K increments
        T   = i_T + T_min;
        // Loop over each velocity
        for (int i_v=0; i_v<=1001; i_v++) {
            // v up in v_esc/1000 increments
            v   = i_v * v_esc/1000;
            A2_Maxwell_cdf[i_T].push_back(Maxwell_flux_cdf(T, v, m));
        }
    }

    return  A2_Maxwell_cdf;
}
//========


double  rand_vel(double T, double xi_1, double xi_2, double T_min,
                 vector< vector<double> > *A2_Maxwell_cdf)
{
    /// Returns a random velocity using the A2_Maxwell_cdf: \
        Finds the 1D subarray for the temperature T, \
        Finds the closest subarray value to the random number xi, \
        Returns a random velocity between the two nearest values in the array.
    double  rand_vel;
    int     v_index,
            T_index = round_to_int(T) - T_min;

    vector<double>::iterator    below_index;

    // Find the first element larger than xi_1
    below_index = lower_bound(((*A2_Maxwell_cdf)[T_index]).begin(),
                              ((*A2_Maxwell_cdf)[T_index]).end(), xi_1);

    v_index     = (below_index - ((*A2_Maxwell_cdf)[T_index]).begin());
    // v goes from 0 to v_esc in 1000 steps
    rand_vel    = (v_index + xi_2) * v_esc/1000;

    return  rand_vel;
}
//========


double  rand_alpha(double xi)
{
    /// Returns a random angle alpha from a random number 0 <= xi < 1 \
        using the inverse cdf for alpha.
    return  asin(xi);
}
//========


double  rand_alpha_flux(double xi)
{
    /// Returns a random angle alpha from a random number 0 <= xi < 1 \
        using the inverse flux cdf for alpha.
    return  0.5 * acos(1 - 2*xi);
}
//========


double  rand_eta(double xi)
{
    /// Returns a random angle eta from a random number 0 <= xi < 1 \
        using the inverse cdf for eta.
    return  2*pi * xi;
}
//========


//============================================================================//
//                              2.2 Migration                                 //
//============================================================================//


double  phi_prime(double phi, double time)
{
    /// Returns the angle phi_ (< 2pi) for a longitude phi at a total time, \
        where the subsolar point is always at phi_ = 0.
    double  phi_;

    phi_    = phi + 2*pi*time/day;
    // Make 0 < phi_ < 2pi
    phi_    = phi_ - 2*pi*floor(phi_/(2*pi));

    return  phi_;
}
//========


double  residence_time(double T, double Q)
{
    /// Returns the residence time from the modified Arrhenius equation.
    return  h/(k*T) * exp(Q * 1e3 / (8.314 * T));
}
//========


//=========================
// Temperature
//=========================
vector< vector<double> >    make_A2_longitude_offset(int seed,
                                                     double sig_offset,
                                                     int num_phi_,
                                                     int num_theta)
{
    /// Make an array of longitude offsets to act as a topography/temperature \
        map for then finding the local temperature in a more physical way.
    double  xi_1, xi_2,
            r, d, n_1,                  // Gaussian RNG
            offset;

    vector< vector<double> >    A2_longitude_offset;

    // Resize array
    A2_longitude_offset.resize(num_phi_);
    for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
        A2_longitude_offset[i_phi_].resize(num_theta);
    }

    // Fill array with Gaussian offsets about zero
    for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
        for (int i_theta=0; i_theta<num_theta; i_theta++) {
            // random normally dist. number using Box-Muller transform
            r   = 0;
            while ((r == 0) || (r > 1)) {
                // -1 < xi < 1
                xi_1    = 2*NR::ran3(seed) - 1;
                xi_2    = 2*NR::ran3(seed) - 1;
                r       = xi_1*xi_1 + xi_2*xi_2;
            }
            d       = sqrt(-2*log(r) / r);
            n_1     = xi_1 * d;

            offset  = sig_offset * n_1;

            A2_longitude_offset[i_phi_][i_theta]    = offset;
        }
    }

    return  A2_longitude_offset;
}
//========


vector< vector<double> >    make_A2_H_temperature(double T_0, double T_1,
                                                  int num_phi_, int num_theta)
{
    /// Makes a 2D array of temperatures for all theta, phi_ [Hurley+2014].
    double  T, phi_, theta,
            theta_s;                    // Solar zenith angle (rad)

    vector< vector<double> >    A2_H_temperature;

    // Resize the array
    A2_H_temperature.resize(num_phi_);
    for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
        A2_H_temperature[i_phi_].resize(num_theta);
    }

    // Set the temperature at each time of day and latitude
    for (int i_phi_=0; i_phi_<num_phi_; i_phi_++) {
        for (int i_theta=0; i_theta<num_theta; i_theta++) {
            phi_    = i_phi_*conv;
            theta   = i_theta*conv;

            // Dayside
            if (is_in_day(phi_)) {
                theta_s = acos(cos(phi_)*cos(pi/2-theta));

                T   = (T_1 + T_0) * pow(cos(theta_s), 0.25);
                if (T < 130) {
                    T   = 130;
                }
            }
            // Nightside
            else {
                T   = 444.738 - 448.937*phi_ + 239.668*pow(phi_,2)
                      - 63.8844*pow(phi_,3) + 8.34064*pow(phi_,4)
                      - 0.423502*pow(phi_,5) + 35*(sin(theta) - 1);
            }

            A2_H_temperature[i_phi_][i_theta]   = T;
        }
    }

    return  A2_H_temperature;
}
//========


//=========================
// Hops
//=========================
double  get_a(double v_0)
{
    /// Returns the semi-major axis a for the elliptical trajectory.
    return  1 / (2/R - v_0*v_0 / (G * M));
}
//========


double  get_e(double a, double alpha)
{
    /// Returns the eccentricity e for the elliptical trajectory.
    double  alpha_  = 0.5*pi - alpha;

    return  sqrt(1 - ((2*a*R-R*R)/(a*a))*sin(alpha_)*sin(alpha_));
}
//========


double  get_max_height(double a, double e)
{
    /// Returns the maximum height reached in the hop.
    return  a*(1+e) - R;
}
//========


double  get_delta(double a, double e)
{
    /// Returns the angle delta used to calculate the end theta coordinate.
    double  cos_delta;

    cos_delta   = (1/e) * (1 - (a * (1 - e*e) / R));

    // Catch nan errors from cos(delta) very slightly above 1 or below -1
    if (cos_delta >= 1) {
        cos_delta   = 1;
    }
    else if (cos_delta <= -1) {
        cos_delta   = -1;
    }

    return  2*acos(cos_delta);
}
//========


double  get_end_theta(double theta_0, double delta, double eta)
{
    /// Returns the theta coordinate after flight.
    return  acos(cos(theta_0)*cos(delta) + sin(theta_0)*sin(delta)*cos(eta));
}
//========


double  get_epsilon(double theta, double theta_0, double delta)
{
    /// Returns the angle epsilon used to calculate the end phi coordinate.
    double  cos_epsilon;

    cos_epsilon = (cos(delta) - cos(theta)*cos(theta_0)) /
                  (sin(theta)*sin(theta_0));

    // Catch nan errors from cos(epsilon) very slightly above 1 or below -1
    if (cos_epsilon >= 1) {
        cos_epsilon = 1;
    }
    else if (cos_epsilon <= -1) {
        cos_epsilon = -1;
    }

    return  acos(cos_epsilon);
}
//========


double  get_end_phi(double phi_0, double epsilon, double eta)
{
    /// Returns the phi coordinate after the hop.
    double phi;

    if (eta > pi) {
        phi = phi_0 + epsilon;
    }
    else {
        phi = phi_0 - epsilon;
    }

    // Keep 0 < phi < 2pi
    if (phi < 0) {
        phi += 2*pi;
    }
    if (phi > 2*pi) {
        phi -= 2*pi;
    }

    return  phi;
}
//========


double  get_flight_time(double a, double e, double delta)
{
    /// Returns the flight time.
    double  T,                          // Orbit period (s)
            E_;                         // Eccentric anomaly (rad)

    T   = sqrt(4*pi*pi*a*a*a/(G*M));

    E_  = pi - 2*atan(sqrt((1-e)/(1+e))*tan((pi-0.5*delta)/2));

    return  (T/pi) * (E_ + e * sin(E_));
}
//========


//=========================
// At altitude
//=========================
double  get_delta_p(double a, double e, double alt)
{
    /// Returns the angle delta_p between the peak and the position at altitude.
    double  cos_delta_p;

    cos_delta_p = 1/e * (1 - (a * (1 - e*e)) / (alt + R));

    // Catch nan errors from cos(delta_p) very slightly above 1 or below -1
    if (cos_delta_p >= 1) {
        cos_delta_p = 1;
    }
    else if (cos_delta_p <= -1) {
        cos_delta_p = -1;
    }

    return  acos(cos_delta_p);
}
//========


double  get_v(double a, double alt)
{
    /// Returns the velocity at altitude.
    return  sqrt(G * M * (2/(R + alt) - 1/a));
}
//========


double  get_alpha(double a, double e, double alt)
{
    /// Returns the angle alpha at altitude.
    double  alpha_;                     // The angle to the vertical (rad)

    alpha_  = asin(sqrt(a*a*(1 - e*e) / (2*a*(R+alt) - (R+alt)*(R+alt))));

    return  90*conv - alpha_;
}
//========


//============================================================================//
//                              2.3 Cold Traps                                //
//============================================================================//


bool    test_capt_at_theta(double theta, double xi, double f_N, double f_S)
{
    /// Returns true if the particle was captured, false if not, depending on a\
        randomly generated number xi.
    double  f;                          // Fractional cold trap area
    bool    is_captured = false;

    if (is_north_pole(theta)) {
        f   = f_N;
    }
    else if (is_south_pole(theta)) {
        f   = f_S;
    }
    else {
        f   = 0;
    }

    if (xi < f) {
        is_captured = true;
    }

    return  is_captured;
}
//========


bool    test_capt_local_traps(double theta, double phi, double theta_trap_N_1,
                              double theta_trap_N_2, double phi_trap_N_1,
                              double phi_trap_N_2, double theta_trap_S_1,
                              double theta_trap_S_2, double phi_trap_S_1,
                              double phi_trap_S_2)
{
    /// Returns true if the particle was captured, false if not, depending only\
        on the location of the particle.
    bool    is_captured = false;

    if ((theta_trap_N_1 < theta) && (theta < theta_trap_N_2) &&
        (phi_trap_N_1 < phi) && (phi < phi_trap_N_2)) {

        is_captured = true;
    }
    else if ((theta_trap_S_1 < theta) && (theta < theta_trap_S_2) &&
             (phi_trap_S_1 < phi) && (phi < phi_trap_S_2)) {
        is_captured = true;
    }

    return  is_captured;
}
//========


bool    test_capt_seasonal(double theta, double xi, double time,
                           double f_season_N, double f_season_S)
{
    /// Returns true if the particle was captured in a seasonal cold trap.
    double  f;
    bool    is_captured = false;

    if (is_north_pole(theta)) {
        f   = f_season_N * sin(2*pi/year * (time - 0.5*year));
    }
    else if (is_south_pole(theta)) {
        f   = f_season_S * sin(2*pi/year * time);
    }
    else {
        f   = 0;
    }

    if (xi < f) {
        is_captured = true;
    }

    return  is_captured;
}
//========


double  get_time_seasonal_rel(double theta, double time, double xi)
{
    /// Returns the time that the particle is released from the seasonal trap.
    int     year_N;                     // Years since start of the simulation
    double  t_trap;                     // Time of capture in seasonal trap (s)

    year_N  = floor(time/year);

    // North
    if (is_north_pole(theta)) {
        // Time since the start of the year at which this trap first appeared
        t_trap  = year/(2*pi) * asin(xi * sin((time - 0.5*year) * 2*pi/year));

        // Release when this same trap disappears
        time    = (year_N + 1)*year - t_trap;
    }
    // South
    else if (is_south_pole(theta)) {
        // Time since the start of the year at which this trap first appeared
        t_trap  = year/(2*pi) * asin(xi * sin(time * 2*pi/year));

        // Release when this same trap disappears
        time    = (year_N + 0.5)*year - t_trap;
    }

    return  time;
}
//========


//============================================================================//
//                              2.4 Misc                                      //
//============================================================================//


int     round_to_int(double x)
{
    /// Returns the integer closest to x.
    return  (x > 0.0) ? (x + 0.5) : (x - 0.5);
}
//========


int     get_index(double val, double max_val, int num_val)
{
    /// Returns the array index for this value.
    return  floor(val/max_val * num_val);
}
//========


double  get_flux(double v, double alpha)
{
    /// Returns the density contribution (``flux'') of the particle. So, weight\
        by the vertical component of the velocity through the surface, i.e. the\
        time spent near the surface.
    return  1 / (v * sin(alpha));
}
//========


//=========================
// Check location
//=========================
bool    is_north_pole(double theta)
{
    /// Returns true if the particle is in the north polar 15 degrees.
    bool    north_pole  = false;

    if ((0*conv < theta) && (theta < 15*conv)) {
        north_pole  = true;
    }

    return  north_pole;
}
//========


bool    is_south_pole(double theta)
{
    /// Returns true if the particle is in the south polar 15 degrees.
    bool    south_pole  = false;

    if ((165*conv < theta) && (theta < 180*conv)) {
        south_pole  = true;
    }

    return  south_pole;
}
//========


bool    is_in_day(double phi_)
{
    /// Returns true if the particle is on the dayside.
    bool    in_day  = false;

    if (not ((90*conv < phi_) && (phi_ < 270*conv))) {
        in_day  = true;
    }

    return  in_day;
}
//========


bool    is_recording_time(double time, int N_bounce, const double start_time,
                          const int start_N_bounce)
{
    /// Returns true if the time and N_bounce are beyond the start values.
    bool    recording_time  = false;

    if ((start_time < time) && (start_N_bounce < N_bounce)) {
        recording_time  = true;
    }

    return  recording_time;
}
//========


bool    is_in_mare(double theta, double phi, double theta_m_1, double theta_m_2,
                   double phi_m_1, double phi_m_2)
{
    /// Returns true if the particle is within the set mare region.
    bool    in_mare = false;
    if (((theta_m_1 < theta) && (theta < theta_m_2)) &&
        ((phi_m_1 < phi) && (phi < phi_m_2))) {

        in_mare = true;
    }

    return  in_mare;
}
//========


bool    is_LADEE_lat(double theta)
{
    /// Returns true if the particle is within +/- 30 degrees of the equator.
    bool    LADEE_lat   = false;

    if ((60*conv < theta) && (theta < 120*conv)) {
        LADEE_lat   = true;
    }

    return  LADEE_lat;
}
//========


//=========================
// Bin Vector Values
//=========================
vector<int> bin_values(vector<double> values, int num_bins, double max_value)
{
    /// Sorts values <= max_value into num_bins bins, \
        returns a vector of the bins.
    int index;

    vector<int> bins;

    bins.resize(num_bins);

    for (int i=0; i<values.size(); i++) {
        if (values[i] <= max_value) {
            index   = floor(values[i]/max_value * num_bins);
            // Put in final bin if value = max_val
            if (index == num_bins) {
                index   = num_bins - 1;
            }
            bins[index] ++;
        }
    }

    return  bins;
}
//========


//============================================================================//
//                              2.5 Admin etc.                                //
//============================================================================//


//=========================
// Check Input Options
//=========================
void    check_arguments(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)
{
    /// Check if any string arguments are invalid, print an error message if so.

    if ((source != "iso_source") && (source != "solar_source") &&
        (source != "local_source") && (source != "cont_iso_source") &&
        (source != "cont_local_source") && (source != "K_source") &&
        (source != "transient_source") && (source != "eqm_source")) {

        cout << "Warning: invalid source!" << endl;
    }

    if ((emerge_dist != "M_B") && (emerge_dist != "M_B_flux")) {
        cout << "Warning: invalid emerge_dist!" << endl;
    }

    if ((losses != "normal") && (losses != "time") &&
        (losses != "normal&time")) {

        cout << "Warning: invalid losses!" << endl;
    }

    if ((T_dist != "H_T_simple") && (T_dist != "H_T_offsets")) {
        cout << "Warning: invalid T_dist!" << endl;
    }

    if ((cold_traps != "prob_at_theta") && (cold_traps != "local_traps") &&
        (cold_traps != "seasonal")) {

        cout << "Warning: invalid cold_traps!" << endl;
    }

    if ((record_dist != "no_dist") && (record_dist != "density") &&
        (record_dist != "densities") && (record_dist != "all_densities") &&
        (record_dist != "at_altitude") && (record_dist != "density_stuck")) {

        cout << "Warning: invalid record_dist!" << endl;
    }

    if ((surf_intn != "Arrhenius") && (surf_intn != "zero") &&
        (surf_intn != "local_Q") && (surf_intn != "squirrelling")) {

        cout << "Warning: invalid surf_intn!" << endl;
    }

    if ((misc != "no_misc") && (misc != "local_Q") &&
        (misc != "mare_T_offset") && (misc != "solar_flare")) {

        cout << "Warning: invalid misc!" << endl;
    }

    if ((rates != "yes_rates") && (rates != "no_rates") &&
        (rates != "LADEE_lat_rates")) {
        cout << "Warning: invalid rates!" << endl;
    }

    if ((backup_results != "yes_backup") && (backup_results != "no_backup")) {
        cout << "Warning: invalid backup_results!" << endl;
    }
}
//========


//=========================
// Printing Information
//=========================
string  get_date_time(char date_y_n, char time_y_n)
{
    /// Returns the current date and/or time.
    time_t  rawtime;
    struct  tm * timeinfo;
    char    buffer[80];

    time(&rawtime);
    timeinfo    = localtime(&rawtime);

    // From arguments, print date and time, just date, or just time
    if ((date_y_n == 'y') && (time_y_n == 'y')) {
        strftime(buffer, 80, "%d,%m,%Y_%H-%M-%S", timeinfo);
    }
    else if ((date_y_n == 'y') && (time_y_n == 'n')) {
        strftime(buffer, 80, "%d/%m/%Y", timeinfo);
    }
    else if ((date_y_n == 'n') && (time_y_n == 'y')) {
        strftime(buffer, 80, "%H:%M:%S", timeinfo);
    }
    // Otherwise, default to both and print error message
    else {
        cout << "Enter 'y' or 'n' to specify date and/or time!" << endl;
        strftime(buffer, 80, "%d/%m/%Y %H:%M:%S", timeinfo);
    }

    string  str(buffer);

    return  str;
}
//========


void    print_date_time(char date_y_n, char time_y_n)
{
    /// Prints the current date and/or time.
    time_t  rawtime;
    struct  tm * timeinfo;
    char    buffer[80];

    time(&rawtime);
    timeinfo    = localtime(&rawtime);

    // From arguments, print date and time, just date, or just time
    if ((date_y_n == 'y') && (time_y_n == 'y')) {
        strftime(buffer, 80, "%d/%m/%Y %H:%M:%S", timeinfo);
    }
    else if ((date_y_n == 'y') && (time_y_n == 'n')) {
        strftime(buffer, 80, "%d/%m/%Y", timeinfo);
    }
    else if ((date_y_n == 'n') && (time_y_n == 'y')) {
        strftime(buffer, 80, "%H:%M:%S", timeinfo);
    }
    // Default to both, print error message
    else {
        cout << "Enter 'y' or 'n' to specify date and/or time!" << endl;
        strftime(buffer, 80, "%d/%m/%Y %H:%M:%S", timeinfo);
    }

    string  str(buffer);

    cout << str << endl;
}
//========


void    print_progress_MPI(int i, int num_loops, int percent, clock_t t_start,
                           int world_rank, int world_size)
{
    /// Prints the percentage completion, the time elapsed, \
        and the estimated time to completion of a loop after the first loop \
        and after each "percent" completion, for a specific MPI processor.
    // No i=0 to avoid division by 0
    if (((i%(num_loops*percent/100) == 0) && (i != 0)) || (i == 1)) {
        clock_t     t_current = clock();
        double      t_elapsed, t_left;

        t_elapsed   = (double) (t_current - t_start) / CLOCKS_PER_SEC;
        t_left      = (double) t_elapsed*num_loops/i - t_elapsed;

        // Print
        cout << "[" << world_rank << "]: Run " << i << " of " << num_loops
             << ", " << (double) 100*i/num_loops << "% complete. ";

        // Print time elapsed in seconds, or minutes or hours if long enough
        if (60*60 < t_elapsed)  {
            cout << round_to_int(t_elapsed/(60*60)) << ":"
                 << round_to_int(fmod(t_elapsed/60, 60)) << ":"
                 << round_to_int(fmod(t_elapsed, 60)) << " elapsed, end in ";
        }
        else if ((60 < t_elapsed) && (t_elapsed < 60*60)) {
            cout << round_to_int(t_elapsed/60) << ":"
                 << round_to_int(fmod(t_elapsed, 60)) << " elapsed, end in ";
        }
        else {
            cout << round_to_int(fmod(t_elapsed, 60)) << "s elapsed, end in ";
        }

        // Print time left in seconds, or minutes or hours if long enough
        if (60*60 < t_left)  {
            cout << round_to_int(t_left/(60*60)) << ":"
                 << round_to_int(fmod(t_left/60, 60)) << ":"
                 << round_to_int(fmod(t_left, 60)) <<  endl;
        }
        else if ((60 < t_left) && (t_left < 60*60)) {
            cout << round_to_int(t_left/60) << ":"
                 << round_to_int(fmod(t_left, 60)) << endl;
        }
        else {
            cout << round_to_int(fmod(t_left, 60)) << "s." << endl;
        }
    }
}
//========




















