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

    Main program.
**/


/*
    0. Includes etc.
    1. Constants
    2. main()
*/


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


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

#include <iostream>
#include <sstream>
#include <ctime>
#include <cmath>
#include <limits>

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. main()                                     //
////////////////////////////////////////////////////////////////////////////////


int main(int argc, char* argv[])
{
    //========
    // Job ID (and index)
    //========
    // Read system argument
    istringstream   ss(argv[1]);
    int job_ID;
    if (!(ss >> job_ID)) {
        cerr << "Invalid job ID " << argv[1] << '\n';
    }

    // Multi job array index (Comment out if not needed; use the placeholder!)
    // Read system argument
    istringstream ss2(argv[2]);
    int job_index;
    if (!(ss2 >> job_index)) {
        cerr << "Invalid job index " << argv[2] << '\n';
    }
    // Placeholder (Must comment out/in if multi job array index is/isn't used!)
//    int job_index   = 0;


    //========
    // Set main simulation parameters
    //========
    // Runs etc.
    int     num_run         = 60*10,        // Total number of simulation runs
            num_particle    = 1e2,          // Number of particles per run
            // Print progress and estimate completion time
            percent         = 20;
            // Catch error for small num_run
            if (num_run*percent/100 < 1) {
                percent     = 100/num_run;
            }

    double  dt              = 5,            // Surface interaction timestep (s)
            // Start and end times (s) and/or number of particle bounces for \
               recording results in the simulation. Set to 0 to ignore.
            start_time      = 200 * day,
            end_time        = 300 * day;
    int     start_N_bounce  = 0,
            end_N_bounce    = 0;

    // Particle processes
    double  mass            = 40,           // Particle mass (amu)
            tau             = 3.14e6,       // Photodestruction timescale (s)
            m               = mass*amu,     // Particle mass (kg)
            Q               = 28,           // Desorption energy (kJ mol^-1)
            P_adsorb        = 1,            // Adsorption probability
            P_squirrel      = 5e-3,         // Squirrelling probability
            f_N             = 4e-2,         // Cold trap fractional areas
                              //1e-4,
            f_S             = f_N;

    // Options
    string  source          = //"iso_source",
                              //"solar_source",
                              //"local_source",
                              //"cont_iso_source",
                              "cont_local_source",
                              //"K_source",
                              //"transient_source",
                              //"eqm_source",

            emerge_dist     = //"M_B",
                              "M_B_flux",

            losses          = //"normal",
                              //"time",
                              "normal&time",

            T_dist          = //"H_T_simple",
                              "H_T_offsets",

            cold_traps      = //"prob_at_theta",
                              //"local_traps",
                              "seasonal",

            record_dist     = //"no_dist",
                              //"density",
                              "densities",
                              //"all_densities",
                              //"at_altitude",
                              //"density_stuck",

            surf_intn       = //"zero",
                              //"Arrhenius",
                              "squirrelling",

            misc            = "no_misc",
                              //"local_Q",
                              //"mare_T_offset",
                              //"solar_flare",

            rates           = "yes_rates",
                              //"no_rates",
                              //"LADEE_lat_rates",

            backup_results  = "yes_backup";
                              //"no_backup";


    //========
    // Run the simulation!
    //========
    run_simulation(num_run, num_particle, percent,
                   dt, start_time, end_time, start_N_bounce, end_N_bounce,
                   mass, tau, m, Q, P_adsorb, P_squirrel, f_N, f_S,
                   source, emerge_dist, losses, T_dist, cold_traps,
                   record_dist, surf_intn, misc, rates, backup_results,
                   job_ID, job_index);


    return  0;
}





















