#include "nuc_eos.h" 
#include "../con2prim.h" 

// constants (required by EOS_Omni)
const double kBerg = 1.38064852e-16;
const double amu   = 1.6605402e-24;

// conversion factors between cgs and M_Sun = c = G = 1
// (temperature between K and MeV)
// see EOS_Omni/doc/units.py
const double rho_gf = 1.61887093132742e-18;
const double press_gf = 1.80123683248503e-39;
const double eps_gf = 1.11265005605362e-21;
const double temp_gf = 8.617330337217213e-11;

// Inverses of the numbers above, calculated manually instead of by
// the compiler
const double inv_rho_gf = 6.17714470405638e17;
const double inv_press_gf = 5.55174079257738e38;
const double inv_eps_gf = 8.98755178736818e20;
const double inv_temp_gf = 1.1604522060401007e10;

// EOS quantities for ideal gas, required by EOS_Omni
double EOS_const[5];

// con2prim parameters
struct c2p_steer c2p;

extern void eos_omni_helmholtzreadtable_(char * );


void collect_metric(struct metric *h, const double hxx, const double hxy, const double hxz,
        const double hyy, const double hyz, const double hzz){

  // get the metric and inverse metric components
  
  for (int i=0;i<4;i++) {
    h->lo[0][i] = 0.0;
    h->lo[i][0] = 0.0;
    h->up[0][i] = 0.0;
    h->up[i][0] = 0.0;
  }
  
  h->lo[1][1] = hxx;
  h->lo[1][2] = hxy;
  h->lo[1][3] = hxz;
  h->lo[2][2] = hyy;
  h->lo[2][3] = hyz;
  h->lo[3][3] = hzz;
        
  h->lo[2][1] = h->lo[1][2];
  h->lo[3][1] = h->lo[1][3];
  h->lo[3][2] = h->lo[2][3];
  
  h->lo_det = hxx*hyy*hzz + 2.0*hxy*hyz*hxz - hxz*hxz*hyy - hxx*hyz*hyz - hxy*hxy*hzz;
  h->lo_sqrt_det = sqrt(h->lo_det);
  h->up_det = 1.0 / h->lo_det;
          
  h->up[1][1] = (-hyz*hyz + hyy*hzz) * h->up_det;
  h->up[1][2] = (hxz*hyz - hxy*hzz) * h->up_det;
  h->up[2][2] = (-hxz*hxz + hxx*hzz) * h->up_det;
  h->up[1][3] = (-hxz*hyy + hxy*hyz) * h->up_det;
  h->up[2][3] = (hxy*hxz - hxx*hyz) * h->up_det;
  h->up[3][3] = (-hxy*hxy + hxx*hyy) * h->up_det;
  
  h->up[2][1] = h->up[1][2];
  h->up[3][1] = h->up[1][3];
  h->up[3][2] = h->up[2][3];
}


int main(int argc, char *argv[]){
  
  if (argc != 7) {
    printf("size of argv= %i\n", argc );
    printf("Usage: name_executable Recovery_scheme EOS param1 param2 Ye rhoT_key");
    printf("Incorrect number of arguments. Terminating...\n");
    exit(1);
  }

  // general EOS settings
  c2p.eoskey_polytrope = 1;
  c2p.eos_prec = 1e-10;
  
  // tolerance for convergence in recovery scheme
  c2p.tol_x = 5e-9;
  c2p.tol_x_retry = 1e2*c2p.tol_x; // reduced tolerance for retry if first attempt failed
  
  // maximum number of iterations for c2p scheme
  c2p.max_iterations   = 300;
  c2p.extra_iterations = 0;
  
  // atmosphere settings
  c2p.rho_atmo = 6e-16;    // in M_Sun = c = G = 1
  c2p.rho_atmo_tol = 0.01; // relative tolerance
  c2p.T_atmo = 9e-07;  // in MeV, 8.62e-06 MeV = 1e5K
  c2p.Ye_atmo = 0.1;   // electron fraction
  c2p.retain_B_atmo = false; // keep the magnetic field in the atmo?
  
  // number of primitive and conserved quantities
  c2p.numprims = 14; // rho, vel123, eps, B123, ye, temp, press, W, ent, abar
  c2p.numcons = 10; // D, S123, tau, B123, ye_con, temp
  
  // stricly enforce v^2 < 1
  c2p.enforce_v2 = 1;
  
  // define relative amplitude for random perturbations 
  // applied to primitives before starting recovery process
  const double perturbation=0.05;
  const double perturbation_W=0.05; // only for Lorentz factor
  
  // tolerance defining successful recovery of primitive variables
  // abs((recovered-original)/original)<tol_prim
  double tol_prim = 10.*c2p.tol_x;
  
  // Choose angle d between B and v
  // See Section 4.1 of Cerda-Duran et al. 2008
  // ALIGNED = 0: d = pi/2
  // ALIGNED = 1: d = 0
  const int ALIGNED = 1;
  
  // EOS parameters for ideal gas and Hybrid EOS
  for(int i=0;i<5;i++)
    EOS_const[i] = 0.0;
  EOS_const[2] = 2.5; //gamma2
  EOS_const[3] = 1.5; //gammath
  EOS_const[4] = ((2.0e14)*rho_gf); //rho_nuc

  // set recovery method
  c2p.c2p_method = atoi(argv[1]);
  c2p.c2p_method_backup = 5;
  c2p.use_c2p_method_backup = false;
  if ((c2p.c2p_method >= 8) || (c2p.c2p_method_backup >= 8)){
    printf("Invalid c2p method(s). Terminating...\n");
    exit(1);
  }

  // set EOS
  c2p.eoskey = atoi(argv[2]);
  if (c2p.eoskey<4) { 
    // ideal gas and Hybrid EOS
    c2p.evolve_T = 0;
    c2p.evolve_Ye = 0;
    EOS_const[0] = 4./3.; 
  } else {
    c2p.evolve_T = 1;
    c2p.evolve_Ye = 1;
  }

  // set electron fraction
  double ye = strtod(argv[5], NULL);

  // set recovery test type
  // plot rho vs. T, or p_mag/p vs. W-1
  double logWminus1 = 0.0;
  double W = 0.0;
  double logPmagP = 0.0;
  double rho_cgs = 0.0;
  double temp_mev = 0.0;
  bool rhoT_key = atoi(argv[6]);
  if (rhoT_key){
    logWminus1 = strtod(argv[3], NULL);
    W = 1.0 + pow(10,logWminus1);
    logPmagP = strtod(argv[4], NULL);
    printf("RS= %i, EOS= %i, logWminus1= %e, logPmagP= %e, ye= %e\n",c2p.c2p_method, c2p.eoskey, logWminus1, logPmagP, ye);
  } else {
    rho_cgs = strtod(argv[3], NULL);
    temp_mev = strtod(argv[4], NULL);
    printf("RS= %i, EOS= %i, rho= %e, temp= %e, ye= %e\n",c2p.c2p_method, c2p.eoskey, rho_cgs, temp_mev, ye);
  }

  // choose size of code test
  // sqrt(npoints) points in each variable (rho,T or W-1,p_mag/p)
  int npoints;
#if DEBUG
  npoints = pow(2,4);
#else
  npoints = pow(2,12);
#endif

  // choose range for variables
  double rho_min, rho_max, temp_min, temp_max;
  // set rho in g/cm^3
  // vary: 1e3 g/cm^3 < rho  < 1e15 g/cm^3
  rho_min = 1.0e3;
  rho_min *= 1.1;  // add 10% to not avoid problems at table boundary
  if (c2p.eoskey==8)
    rho_max = 9.0e11;
  else
    rho_max = 1.0e15;

  // set temp in K
  // vary 1e5*Kelvin < temp < 50 MeV
  if (c2p.eoskey==4) {
    temp_min = 1.16e8;
    temp_min *= 1.1;  // add 10% to not avoid problems at table boundary
  } else {
    temp_min = 1e5;
  }    
  temp_max = 50.0*inv_temp_gf;

  // set Wminus1, PmagP
  double logWminus1_min, logWminus1_max, logPmagP_min, logPmagP_max;
  // vary W-1 from 1e-4 to 1e4
  logWminus1_min = -4.0;
  logWminus1_max = 4.0;
  // vary P_mag/P from 1e-8 to 1e10
  logPmagP_min = -8.0;
  logPmagP_max = 10.0;

  printf("rho_min = %g\n", rho_min);
  printf("rho_max = %g\n", rho_max);
  printf("temp_min = %g\n", temp_min);
  printf("temp_max = %g\n", temp_max);
  printf("logWminus1_min = %g\n", logWminus1_min);
  printf("logWminus1_max = %g\n", logWminus1_max);
  printf("logPmagP_min = %g\n", logPmagP_min);
  printf("logPmagP_max = %g\n", logPmagP_max);

  // READ EOS TABLES
  
  // tabulated nuclear EOS
  if(c2p.eoskey == 4){
    nuc_eos_C_ReadTable("eostable.h5");
    printf("nrho: %d\n",nrho);
    printf("ntemp: %d\n",ntemp);
    printf("nye: %d\n\n",nye);
  }

  // Helmholtz EOS
  if(c2p.eoskey == 8){
    char helmeos_table_name[50];
    strcpy(helmeos_table_name, "helm_table_FLASH.dat");
    eos_omni_helmholtzreadtable_(helmeos_table_name);
  }

  // further EOS settings
  double rhomi,rhoma,tmi,tma,yemi,yema,emi,ema,pmi,pma;
  EOS_Omni_get_eos_limits(c2p.eoskey,&rhomi,&rhoma,&tmi,&tma,&yemi,&yema,&emi,&ema,&pmi,&pma);
  c2p.eos_rho_min = rhomi;
  c2p.eos_rho_max = rhoma;
  c2p.eos_temp_min = tmi;
  c2p.eos_temp_max = tma;
  c2p.eos_eps_min = emi; //-1.518793e-03; //emi;
  c2p.eos_eps_max = ema;
  c2p.eos_press_min = pmi;
  c2p.eos_press_max = pma;
  printf("rma = %g,rmi = %g,temp_ma = %g,temp_mi = %g,eps_ma = %g,eps_mi = %g,pmi=%g,pma=%g\n",rhoma,rhomi,tma,tmi,ema,emi,pmi,pma);

  // Open output file

  FILE *fp;
  char fname[256];
  if (rhoT_key)
    sprintf(fname, "c2p_output_A%i_RS-%i_EOS-%i_logWminus1-%.2g_logPmagP-%.2g_YE-%.1f.dat", ALIGNED, c2p.c2p_method, c2p.eoskey, logWminus1, logPmagP, ye);
  else
    sprintf(fname, "c2p_output_A%i_RS-%i_EOS-%i_RHO-%.2g_T-%.2g_YE-%.1f.dat", ALIGNED, c2p.c2p_method, c2p.eoskey, rho_cgs, temp_mev, ye);
  
  fp = fopen(fname, "w");
  if (fp == NULL) {
    printf("Couldn't open output file for writing.\n");
    exit(1);
  }
  
  // Write setup parameters to output file

  if(c2p.eoskey==1)
    fprintf(fp,"# Polytropic EOS\n");
  else if(c2p.eoskey==2)
    fprintf(fp,"# Ideal gas EOS\n");
  else if(c2p.eoskey==3)
    fprintf(fp,"# Hybrid EOS\n");
  else if(c2p.eoskey==4)
    fprintf(fp,"# Tabulated EOS\n");
  else if(c2p.eoskey==8)
    fprintf(fp,"# Helmholtz EOS\n");  
  else{
    printf("EOS not implemented. Terminating...\n");
    exit(1);
  }

  if(ALIGNED==0)
    fprintf(fp,"# NOT ALIGNED: cos^2(delta) = 0\n");
  else if(ALIGNED==1)
    fprintf(fp,"# ALIGNED: cos^2(delta) = 1\n");
  else {
    printf("ALIGNED = %i parameter not known. Aborting.", ALIGNED);
    exit(1);
  }

  if (rhoT_key) {
    fprintf(fp,"# rhoT_key = 1");
    fprintf(fp,"# logWminus1 = %e\n",logWminus1);
    fprintf(fp,"# logPmagP = %e\n",logPmagP);
  } else {
    fprintf(fp,"# rhoT_key = 0");
    fprintf(fp,"# rho_cgs = %e\n",rho_cgs);
    fprintf(fp,"# temp_mev = %e\n",temp_mev);
  }

  // define the metric
  struct metric g;
  collect_metric(&g, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0);

  if(g.lo_det <= 0.0){
    printf("det(g) <= 0. Con2prim does not make sense at this point.");
  }
  
  // allocate primitive and conservative arrays
  double* prim = malloc(sizeof(double) * c2p.numprims);
  double* prim_all = malloc(sizeof(double) * npoints * c2p.numprims);
  double* original_prim = malloc(sizeof(double) * npoints * c2p.numprims);
  if (prim_all == NULL) {
	  printf("Could not allocate enough much memory for prims...\n");
	  exit(1);
  }
  double* con = malloc(sizeof(double) * c2p.numcons);
  double* con_all= malloc(sizeof(double) * npoints * c2p.numcons);
  if (con_all == NULL) {
	  printf("Could not allocate enough memory for cons...\n");
	  exit(1);
  }
  
  // initialize prims and cons
  for(int i=0;i<(npoints*c2p.numprims);i++){
    prim_all[i]=0.0;
  }
  for(int i=0;i<(c2p.numprims);i++){
    prim[i]=0.0;
  }
  for(int i=0;i<(npoints*c2p.numcons);i++){
    con_all[i]=0.0;
  }
  for(int i=0;i<(c2p.numcons);i++){
    con[i]=0.0;
  }

  //------------------------------------------------------------
  // Generate grid of primitive quantities and compute conserved
  // quantities for subsequent inversion 
  //------------------------------------------------------------

  double n1,n2;
  double v,vx,vy,vz;
  double rhoexp;
  double delrho;
  double exp1;
  double del1;
  double tempexp;
  double deltemp;
  double tempexp_i;
  double rhoexp_i;

  // define range of rho or W-1 depending on rhoT_key
  if (rhoT_key) {
    rhoexp_i = log10(rho_min);
    rhoexp = rhoexp_i;
    delrho = (log10(rho_max)-rhoexp_i)/(sqrt(1.0*npoints)-1.0);
  } else {
    exp1 = logWminus1_min;
    del1 = (logWminus1_max-logWminus1_min)/(sqrt(1.0*npoints)-1.0);
  }

  printf("\n---------------------------------------------------\n");
  printf("Generating conservatives: prim2con\n");
  printf("---------------------------------------------------\n\n");
  
  for(n1=0;n1<sqrt(1.0*npoints);n1++){
    
    double exp2;
    double del2;

    if (rhoT_key) {
      tempexp_i = log10(temp_min*temp_gf);
      tempexp = tempexp_i;
      deltemp = (log10(temp_max*temp_gf)-tempexp_i)/(sqrt(1.0*npoints)-1.0);
    } else {
      exp2 = logPmagP_min;
      del2 = (logPmagP_max-logPmagP_min)/(sqrt(1.0*npoints)-1.0);
    }

    for(n2=0;n2<sqrt(1.0*npoints);n2++){
      
      // Translate 2D indexing to 1D index
      int i = n1*sqrt(1.0*npoints)+n2;
	
	    printf("**************************************************\n"); 
      printf("Creating data set No.: %d\n", i+1); 
      printf("n1,n2 = %g,%g\n", n1,n2);
      if (rhoT_key) {
        rho_cgs  = pow(10, rhoexp);
        temp_mev = pow(10,tempexp);
      } else {
        // (W-1) = 10^(exp1)
        W = 1.0+pow(10.0,exp1);
      }

      // Velocity magnitude
      v = sqrt(1.0-1.0/(W*W));
      // Randomly choose vx, vy, vz
      vx = v*((double)rand())/((double)RAND_MAX);
      vy = sqrt(v*v - vx*vx)*((double)rand())/((double)RAND_MAX);
      vz = sqrt(v*v - vx*vx - vy*vy);
      // Verify the magnitude
      if(fabs(v*v - (vx*vx+vy*vy+vz*vz))>1.0e-10){
  	    printf("v problem!!!\n");
  	    exit(1);
      }
      
      // Set primitive variables
      prim_all[i*c2p.numprims+RHO] = rho_cgs*rho_gf;
      prim_all[i*c2p.numprims+v1_cov] = vx;
      prim_all[i*c2p.numprims+v2_cov] = vy;
      prim_all[i*c2p.numprims+v3_cov] = vz;
      prim_all[i*c2p.numprims+TEMP] = temp_mev;
      prim_all[i*c2p.numprims+YE] = ye;
      prim_all[i*c2p.numprims+WLORENTZ] = W;
	  
  	  for(int k=0;k<c2p.numprims;k++){
  	    prim[k] = prim_all[i*c2p.numprims+k];
  	  }
 	  
      // compute eps, press
      double epsilon, pressure;
  	  epsilon = eps_from_prim(prim);
  	  prim_all[i*c2p.numprims+EPS] = epsilon; 
  	  prim[EPS] = epsilon;
      pressure = P_from_prim(prim);
  	  prim_all[i*c2p.numprims+PRESS] = pressure;
  	  prim[PRESS] = pressure;
        
      // check P > 0
      if(pressure<0.0){
  	    printf("negative initial pressure: press = %e\n",pressure);
  	    printf("Terminating...\n");
  	    exit(1);
      }
      
      // set up magnetic field
      double B;
      if (rhoT_key) {
        B = sqrt(2.0*pow(10.0,logPmagP)*pressure);
      } else {
        B = sqrt(2.0*pow(10.0,exp2)*pressure);
      }
      double Bxhat;
      double Byhat;
      double Bzhat;
      if(ALIGNED==1){

  	    double rand_temp = ((double)rand())/((double)RAND_MAX);
  	
  	    if(rand_temp<0.5){
  	      Bxhat = vx/v;
  	      Byhat = vy/v;
  	      Bzhat = vz/v;
  	    }
  	    else{
  	      Bxhat = -1.0*vx/v;
  	      Byhat = -1.0*vy/v;
  	      Bzhat = -1.0*vz/v;
  	    }
      }
      else{
  	    Bxhat = ((double)rand())/((double)RAND_MAX);
  	    Byhat = sqrt(1.0 - Bxhat*Bxhat)*((double)rand())/((double)RAND_MAX);
  	    Bzhat = (-Bxhat*vx-Byhat*vy)/vz;
  	    double normB = sqrt(Bxhat*Bxhat + Byhat*Byhat + Bzhat*Bzhat);
  	    Bxhat = Bxhat/normB;
  	    Byhat = Byhat/normB;
  	    Bzhat = Bzhat/normB;
  	
  	    // Verify orthogonality
  	    if(fabs((Bxhat*vx+Byhat*vy+Bzhat*vz))>1.0e-10){
  	      printf("Orthogonality problem!\n");
  	      printf("Terminating...\n");
  	      exit(1);
  	    }
      }
      prim_all[i*c2p.numprims+B1_con] = Bxhat*B;
      prim_all[i*c2p.numprims+B2_con] = Byhat*B;
      prim_all[i*c2p.numprims+B3_con] = Bzhat*B;
      // Verify the magnitude
      if(fabs(1.0 -(Bxhat*Bxhat+Byhat*Byhat+Bzhat*Bzhat))>1.0e-10){
        printf("B problem!!!\n");
        printf("%e\n",Bxhat*Bxhat+Byhat*Byhat+Bzhat*Bzhat);
        printf("Terminating...\n");
        exit(1);
      }
      
      // Save original primitive variables
      // and set primitive variables for present prim2con conversion
      for(int k=0;k<c2p.numprims;k++){
  	    original_prim[i*c2p.numprims+k] = prim_all[i*c2p.numprims+k];
  	    prim[k] = prim_all[i*c2p.numprims+k];
        printf("Save prims:\n");
        printf("%e\t%e\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k]);
      }
      
      // Set up conservatives by converting prim to con
      prim2con_MHD_(prim, con, g.up, g.lo);
      // update con_all
      for(int k=0;k<c2p.numcons;k++){
        con_all[i*c2p.numcons+k] = con[k];
      }

#if DEBUG
      // print primitive and conservative data sets
      printf("\nprim2con Results (i=%d):\nprims\t\tcons\n",i);
      for(int k=0;k<c2p.numprims;k++){
        if (k>=c2p.numcons) {
          printf("%e\t%e\n",prim[k],0.0);
        } else{
          printf("%e\t%e\n",prim[k],con[k]);
        }  
      }

      for(int k=0;k<c2p.numcons;k++){
        printf("final2 p2c:\n");
        printf("con[%i] \t con_all[%i]\n", k, i*c2p.numcons+k);
        printf("%g \t %g\n", con[k], con_all[i*c2p.numcons+k]);
      }
#endif
      if (rhoT_key) {
        tempexp = tempexp + deltemp;
      } else {
        exp2 = exp2 + del2;
      }

    }
    if (rhoT_key) {
      rhoexp = rhoexp + delrho;
    } else {
      exp1 = exp1 + del1;
    }

  }
  
  //------------------------------------------------------------
  // Conservative to Primitive Recovery
  //------------------------------------------------------------

  // Array to store number of iterations required for con2prim
  // at each point
  int count[npoints];  // number of iterations
  int nEOScalls[npoints];  // number of EOS calls
  double error[npoints];   // relative accuracy
  bool retry[npoints]; // record whether c2p has been carried out with reduced tolerance
  int failed[npoints]; // record whether c2p has failed
  int c2p_keyerr[npoints]; // record c2p error code
  
  // Reset primitive variables/setting up initial guesses
  for(int i=0;i<npoints;i++){
    for(int k=0;k<c2p.numprims;k++){
      int randnum = ((int)rand() % 2); // 0 or 1
      randnum = randnum*2; // 0 or 2
      randnum = randnum - 1; // -1 or 1
      if ((k == RHO) && (c2p.c2p_method == 3)){
        prim_all[i*c2p.numprims+k] *= (1.0 + ((double)randnum)*perturbation_W);
      } else if (k == WLORENTZ){
        // perturb W-1
        prim_all[i*c2p.numprims+k] = 1.0 + (prim_all[i*c2p.numprims+k] - 1.0) * (1.0 + ((double)randnum)*perturbation_W);
      } else {
        // perturb all other primitives
        prim_all[i*c2p.numprims+k] *= (1.0 + ((double)randnum)*perturbation);
      }
    }

    // checks on initial guesses
    // ensure they are physical
    
    // density
    if (prim_all[i*c2p.numprims+RHO] < c2p.eos_rho_min){
      prim_all[i*c2p.numprims+RHO] = (1.0+perturbation)*c2p.eos_rho_min;
    }
    if (prim_all[i*c2p.numprims+RHO] > c2p.eos_rho_max){
      prim_all[i*c2p.numprims+RHO] = (1.0-perturbation)*c2p.eos_rho_max;
    }
    
    // velocities
    // make sure that initial guesses satisfy v<c
    vx = prim_all[i*c2p.numprims+v1_cov];
    vy = prim_all[i*c2p.numprims+v2_cov];
    vz = prim_all[i*c2p.numprims+v3_cov];
    if(sqrt(vx*vx+vy*vy+vz*vz) >= 1.0){
	    printf("v problem with initial guesses...!\n");
	    printf("...reset to original values!\n");
	    prim_all[i*c2p.numprims+v1_cov] = original_prim[i*c2p.numprims+v1_cov];
	    prim_all[i*c2p.numprims+v2_cov] = original_prim[i*c2p.numprims+v2_cov];
	    prim_all[i*c2p.numprims+v3_cov] = original_prim[i*c2p.numprims+v3_cov];
    
      vx = prim_all[i*c2p.numprims+v1_cov];
      vy = prim_all[i*c2p.numprims+v2_cov];
      vz = prim_all[i*c2p.numprims+v3_cov];
      if(sqrt(vx*vx+vy*vy+vz*vz) >= 1.0){
	      printf("Reset of v did not work... something went wrong. Terminating...\n");
        exit(1);
      }
    }
    
    // specific internal energy
    if (prim_all[i*c2p.numprims+EPS] < c2p.eos_eps_min){
      prim_all[i*c2p.numprims+EPS] = fabs((1.0+perturbation)*c2p.eos_eps_min - c2p.eos_eps_min) + c2p.eos_eps_min;
    }
    if (prim_all[i*c2p.numprims+EPS] > c2p.eos_eps_max){
      prim_all[i*c2p.numprims+EPS] = (1.0-perturbation)*c2p.eos_eps_max;
    }
  }

  // Start the con2prim process
  printf("\n---------------------------------------------------\n");
  printf("Starting con2prim\n");
  printf("---------------------------------------------------\n\n");

  for(int i=0;i<npoints;i++){
    printf("**************************************************\n");
    printf("Starting inversion No. %d\n", i+1);
    
    // set primitives (initial guesses) and conservatives
    for(int k=0;k<c2p.numprims;k++){
	    prim[k] = prim_all[i*c2p.numprims+k];
    }
    for(int k=0;k<c2p.numcons;k++){
      con[k] = con_all[i*c2p.numcons+k];
    }

    // recover Ye from conservatives
    prim[YE]=con[YE]/con[D];
    
#if DEBUG    
    // check initial guesses
    printf("primitives and conservatives:\n");
    printf("prim guesses\tprims\t\tcons\n");
    for(int k=0;k<c2p.numprims;k++)
    {
      if (k< c2p.numcons)
      {
        printf("%e\t%e\t%e\n",prim[k],original_prim[i*c2p.numprims+k],con_all[i*c2p.numcons+k]);
      } else {
        printf("%e\t%e\n",prim[k],original_prim[i*c2p.numprims+k]);
      }
    }
#endif

    // con2prim process
    struct c2p_report report;
    bool excise = false;
    bool c2p_grace = false;
          
    failed[i] = 0;
    c2p_keyerr[i] = 0;
    retry[i] = false;
    // call con2prim process
    con2prim_MHD_(prim,con,g.up,g.lo,excise,c2p_grace,&report);
    count[i] = report.count;
    retry[i] = report.retry;
    nEOScalls[i] = report.nEOScalls;
    
    if (report.failed) {
      failed[i] = 1;
      c2p_keyerr[i] = report.c2p_keyerr;
      printf("FAILED: con2prim reconstruction failed for the following reason:\n");
      printf("Report: %s\n", report.err_msg);
    } else {
      printf("SUCCESS: con2prim converged in %i iterations!\n",count[i]);
    }
    
    // update primitives
    for(int k=0;k<c2p.numprims;k++){
	    prim_all[i*c2p.numprims+k] = prim[k];
    }
  }
  
  // free storage of EOS tables
  if(c2p.eoskey==4){
    free(alltables);
    free(epstable);
    free(presstable);
    free(logrho);
    free(logtemp);
    free(yes);
  }


  //------------------------------------------------------------
  // Analysis
  //------------------------------------------------------------

  // check for convergence and successful recovery of primitives
  // count number of iterations for successful and failed recoveries
  double nofailaverage = 0.0;
  double average = 0.0;
  int nofailtotal = 0;
  int nofailconv = 0;
  int converged;
  int recovery;
  double tol_prim_loc;

  
  printf("\n---------------------------------------------------\n");
  printf("Analysis of con2prim\n");
  printf("---------------------------------------------------\n\n");
  
  for(int i=0;i<npoints;i++){
    // For each point, check for discrepancy in all primitive variables
    // check that appropriate tolerance is used
    if (retry[i]){
      tol_prim_loc = tol_prim*c2p.tol_x_retry/c2p.tol_x;
    }else{
      tol_prim_loc = tol_prim;
    }
    // check for discrepancy  
    recovery = 1;
    error[i] = 0.0;
    for(int k=0;k<c2p.numprims;k++){
      if((!(c2p.eoskey==1 && k == EPS))&& ( (k==RHO) || (k==v1_cov) || (k==v2_cov) || (k==v3_cov) || (k==EPS) ))
      { 
        //tolerance enforced on rho, vel, eps

  	    if((fabs((original_prim[i*c2p.numprims+k] -
  	      prim_all[i*c2p.numprims+k])/original_prim[i*c2p.numprims+k])>tol_prim_loc) ||
  	      (prim_all[i*c2p.numprims+k]!=prim_all[i*c2p.numprims+k])){
  	      // If there is a discrepancy
#if DEBUG          
  	      printf("i,k, discrepancy: %i, %i, %e\n", i,k,fabs((original_prim[i*c2p.numprims+k] -
  	      prim_all[i*c2p.numprims+k])/original_prim[i*c2p.numprims+k]));
          printf("%g\t%g\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k]);
#endif
  	      recovery = 0;
  	    }
        error[i] += fabs((original_prim[i*c2p.numprims+k]-prim_all[i*c2p.numprims+k]) / original_prim[i*c2p.numprims+k]);
      }
    }

    if (!(c2p.eoskey==1)) {
      error[i] = 0.2 * error[i];
    } else {
      error[i] = 0.25 * error[i];
    }  
    // count the total number of iterations
    average += count[i];

    // check if scheme converged and count the corresponding number of iterations
    if (failed[i] != 1) {
      nofailconv++;
    }
    
    // count the successful recoveries and the corresponding number of iterations
    if(recovery==1){
      nofailtotal++;
      nofailaverage += count[i];
    }

    // check if eps<0 and v>c
    int epsneg = 0;
    if(prim_all[i*c2p.numprims+EPS]<0.0)
      epsneg = 1;
    int gtc = 0;
    if((sqrt(prim_all[i*c2p.numprims+v1_cov]*prim_all[i*c2p.numprims+v1_cov]+prim_all[i*c2p.numprims+v2_cov]*prim_all[i*c2p.numprims+v2_cov]+prim_all[i*c2p.numprims+v3_cov]*prim_all[i*c2p.numprims+v3_cov]) -1.0)>=0.0)
      gtc = 1;
    
    // for each point, store its parameters (log(rho), log(temp)) or (log(W-1), log(Pmag/P))
    // an recovery characteristics in output file
    int n2 = (i%((int)sqrt(1.0*npoints)));
    int n1 = (i-n2)/sqrt(1.0*npoints);
    double del1;
    double del2;
    double exp1;
    double exp2;

    if (rhoT_key)
    {
      tempexp = tempexp_i+n2*deltemp;
      rhoexp = rhoexp_i + n1*delrho;
    } else {
      del1 = (logWminus1_max-logWminus1_min)/(sqrt(1.0*npoints)-1.0);
      del2 = (logPmagP_max-logPmagP_min)/(sqrt(1.0*npoints)-1.0);
      exp1 = logWminus1_min+n1*del1;
      exp2 = logPmagP_min+n2*del2;
    }

    if (rhoT_key) {
      fprintf(fp,"%e %e %d %i %i %d %d %i %e %e %e %i\n",logWminus1,logPmagP,recovery,failed[i],count[i],epsneg,gtc, c2p_keyerr[i], rhoexp, tempexp, error[i], nEOScalls[i]);
    } else {
      fprintf(fp,"%e %e %d %i %i %d %d %i %e %e %e %i\n",exp1,exp2,recovery,failed[i],count[i],epsneg,gtc, c2p_keyerr[i], rhoexp, tempexp, error[i], nEOScalls[i]);
    }

#if DEBUG
    // print recovery results
    printf("**************************************\n");
    printf("Con2prim results, inversion No.: %d\n", i+1);
    if(recovery==0)
    {
      if (rhoT_key) {
        printf("\n*****FAILED recovery at\nrho\ttemp\n%e\t%e\nafter %d iterations\n\n",pow(10,rhoexp), pow(10,tempexp),count[i]);
      } else {
        printf("\n*****FAILED recovery at\nlog(W-1)\tlog(Pmag/P)\n%e\t%e\nafter %d iterations\n\n",exp1,exp2,count[i]);
      }
    }
    else
    { 
      if (rhoT_key) {
      printf("\nSUCCESSFUL recovery at\nrho\ttemp\n%e\t%e\nafter %d iterations\n\n",pow(10,rhoexp), pow(10,tempexp),count[i]);
      } else {
      printf("\nSUCCESSFUL recovery at\nlog(W-1)\tlog(Pmag/P)\n%e\t%e\nafter %d iterations\n\n",exp1,exp2,count[i]);
      }
    }
    printf("recovered prims\toriginal prims\tcons\n");  
    for(int k=0;k<c2p.numprims;k++)
    {
      if (k< c2p.numcons) {
        printf("%e\t%e\t%e\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k],con_all[i*c2p.numcons+k]);
      } else {
        printf("%e\t%e\n",prim_all[i*c2p.numprims+k],original_prim[i*c2p.numprims+k]);
      }
    }
#endif

  }   
    
  // print final report 
  printf("\nFINAL REPORT:\n"); 
  printf("%d/%d recovered (%f percent)\n",nofailtotal,npoints,nofailtotal/(1.0*npoints)*100.0);
  printf("%d/%d converged (%f percent)\n",nofailconv,npoints,nofailconv/(1.0*npoints)*100.0);   
  printf("Average number of iterations (all): %f\n",average/(1.0*npoints));
  printf("Average number of iterations (recovered only): %f\n\n",nofailaverage/(1.0*nofailtotal)); 
  if (rhoT_key==1){
    printf("W=%e, logPmag/P= %e\n",W, logPmagP );
  } else {
    printf("rho=%e, temp_mev= %e\n",rho_cgs, temp_mev );
  }
  // free arrays
  free(prim);
  free(con);
  free(prim_all);
  free(con_all);
  free(original_prim);
  
  return 0;
}