#include "iteration_functions.h"
#include "globals.h"
#include <mpi.h>
#include <stdio.h>

char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method);

int seed;
double tau_c;
double albedobase;


char* polarization_file;
char* diffusion_file;


//int nph;
//int rank;

int kiter;
int DiskPolarizationFlag;
int Albedoflag;

//char* outspec;
//char* polarfile;
char* integral;
//char* diffusion;
char* polarinput;

FILE* dat_spectrum;
FILE* dat_diffusion;
FILE* dat_integral_polar;
FILE* dat_energy_polar;

int main(int argc, char* argv[])
{
  int t;
  int status;
  char* method;
  FlagEleDist = MAXW;

  /*======================================================================*/
  /*Initialize MPI*/
  /*======================================================================*/

  if (argc == 1)
  {
    printf(
        "Program to compute with an iterative method the polarization degree and angular "
        "distribution at the top of a plane-parallel scattering atmosphere\n");

    printf("Synopsis\n");
    printf("start_interation  disktau=disktau seed=[1,2,3] method=ode|st85|mc\n\n");

    printf("Shared parameters\n\n");
    printf("disktau=vaue: slab optical depth\n\n");

    printf("\nSeed photons distributions: three options are possible with parameter seed=N\n");
    printf("1: photons at the center of the slab\n");
    printf("2: photons uniform distribution\n");
    printf("3: photons at the base (valid only for ode and mc methods)\n\n");

    printf("Parameters for iteration algorithms\n");

    printf("method=ode  : solution of the ODE system for Il and Ir\n");
    printf(
        "method=st85 : use ST85 algorithm (valid only for seed photon at the center or uniformly "
        "distributed)\n");

    printf("kiter=N: number of iterations\n\n");

    printf("Parameters for MC simulation\n");
    printf("method=mc\n");
    printf("nph=N       : number of simulated photons\n");
    printf("kte=value   : electron temperature (default is 1 keV)\n");
    printf("For seed photons at the bottom of the slab select if they are unpolarized\n");
    printf(
        "with option polardisk=n or the polarization and angular distribution is read from file\n");
    printf("with polardisk=file\n");

    return (1);
  }

  char* tau_char;
  char* seed_char;
  char* albedo_char;

  for (t = 1; t < argc; t++)
  {
    if (strstr(argv[t], "disktau=") != NULL)
    {
      disktau = numpar(argv[t]);
      tau_char = stringpar(argv[t]);
    }

    if (strstr(argv[t], "kiter=") != NULL)
      kiter = (int)numpar(argv[t]);

    if (strstr(argv[t], "seed=") != NULL)
    {
      seed = (int)numpar(argv[t]);
      seed_char = stringpar(argv[t]);
    }

    if (strstr(argv[t], "tau_c=") != NULL)
      tau_c = numpar(argv[t]);

    if (strstr(argv[t], "method=") != NULL)
    {
      method = stringpar(argv[t]);
    }

    if (strstr(argv[t], "nph=") != NULL)
    {
      nph = (int)numpar(argv[t]);
    }

    if (strstr(argv[t], "polarfile=") != NULL)
    {
      polarfile = stringpar(argv[t]);
    }

    if (strstr(argv[t], "diffusion=") != NULL)
    {
      diffusion = stringpar(argv[t]);
    }

    if (strstr(argv[t], "outspec=") != NULL)
    {
      outspec = stringpar(argv[t]);
    }

    if (strstr(argv[t], "integral=") != NULL)
    {
      integral = stringpar(argv[t]);
    }

    if (strstr(argv[t], "polardisk=") != NULL)
    {
      polardisk = stringpar(argv[t]);
    }

    if (strstr(argv[t], "polarinput=") != NULL)
    {
      polarinput = stringpar(argv[t]);
    }

    if (strstr(argv[t], "kte=") != NULL)
    {
      kte = numpar(argv[t]);
    }

    if (strstr(argv[t], "ktbb=") != NULL)
    {
      ktbb = numpar(argv[t]);
    }

    if (strstr(argv[t], "albedobase=") != NULL)
    {
      albedobase = numpar(argv[t]);
      albedo_char = stringpar(argv[t]);
      Albedoflag = 1;
    }
  }

  /*=========================================================*/
  /*Default values*/
  /*=========================================================*/

  if (method == NULL || method == 0x0)
  {
    printf("Please select the method to compute results with method=ode|st85|mc\n");
    return 1;
  }

  if (!kte)
  {
    kte = 1;
  }

  if (!ktbb)
  {
    ktbb = 1;
  }

  /*=========================================================*/

  if (!disktau || disktau == 0)
  {
    printf("Please select the slab optical depth with parameter disktau=value\n");
    return (0);
  }

  if (!seed)
  {
    printf("Please select the seed photon distribution with parameter seed=1,2,3\n\n");
    return (1);
  }

  if (strcmp(method, "ode") == 0)
  {
    diffusion_file = set_filename("diffusion_ode_tau", tau_char, seed_char, albedo_char, "ode");
    polarization_file = set_filename("intensity_ode_tau", tau_char, seed_char, albedo_char, "ode");

    if (kiter == 0)
    {
      printf("Please select the number of iterations with parameter kiter=kiter\n");
      exit(1);
    }

    optimize_taustep(kiter, seed);
  }
  else if (strcmp(method, "st85") == 0)
  {
    polarization_file = set_filename("intensity_st85_tau", tau_char, seed_char, albedo_char, "st85");
    diffusion_file = set_filename("diffusion_st85_tau", tau_char, seed_char, albedo_char, "st85");

    if (kiter == 0)
    {
      printf("Please select the number of iterations with parameter kiter=kiter\n");
      exit(1);
    }

    slab_polarization(kiter, seed);
  }
  else if (strcmp(method, "mc") == 0)
  {
    MPI_Init(&argc, &argv);
    init_mpi();

    if (polardisk == 0x0 && seed == 3 && rank == 0)
    {
      printf("Set polardisk=n or polardisk=file\n");
      exit(1);
    }

    if (seed == 1 || seed == 2)
    {
      polardisk = "n";
    }

    if (Albedoflag != 1)
    {
      printf("Please select the albedo at the base of the slab with albedobase=value\n");
      exit(1);
    }

    if (albedobase < 0 || albedobase > 1)
    {
      printf("Error: albedo must be in the range 0-1\n");
      exit(1);
    }

    if (nph == 0)
    {
      printf("Please select the number of photons for MC simulation with parametere nph=nph\n");
      exit(1);
    }

    if (strcmp(polardisk, "n") == 0)
    {
      DiskPolarizationFlag = FALSE;
    }
    else if (strcmp(polardisk, "file") == 0)
    {
      DiskPolarizationFlag = FROMFILE;

      if (polarinput == 0x0)
      {
        printf("Select the pre-tabulated reference file with parameter polarinput=filename\n");
        exit(1);
      }
    }
    else
    {
      printf("Unknonw option for parameter polardisk\n");
      exit(1);
    }

    diffusion = set_filename("diffusion_mc_tau", tau_char, seed_char, albedo_char, "mc");
    outspec = set_filename("spectrum_mc_tau", tau_char, seed_char, albedo_char, "mc");
    integral = set_filename("integralpolar_mc_tau", tau_char, seed_char, albedo_char, "mc");
    polarfile = set_filename("energypolar_mc_tau", tau_char, seed_char, albedo_char, "mc");

    //status = slab_mc(nph, seed);

    simul_sync();
    
    
    // check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);
  }

  else
  {
    printf(
        "\nPlease select method=ode (solve ODE), method=st85 (ST85) or method=mc (Montecarlo)\n");
    return (1);
  }

  if (strcmp(method, "mc") == 0) MPI_Finalize();

  return (0);
}

char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method)
{
  // root="intensity_ode";

  char* filename;
  char* tmp_a;
  char* tmp_b;
  char* tmp_c;
  char* tmp_d;

  tmp_a = stringconcat(root, tau_char);

  // printf("tmp_a = %s\n", tmp_a);

  tmp_b = stringconcat("_seed", seed_char);

  // printf("tmp_b = %s\n", tmp_b);

  tmp_c = stringconcat(tmp_a, tmp_b);

  // printf("tmp_c = %s\n", tmp_c);

  if (strcmp(method, "mc") == 0)
  {
    tmp_d = stringconcat(tmp_c, stringconcat("_A", albedo));

    filename = stringconcat(tmp_d, ".qdp");
  }
  else
  {
    filename = stringconcat(tmp_c, ".qdp");
  }

  // printf("Final name %s\n\n", filename);

  return filename;
}

/*==============================================================================*/