/*@@
   @file      EOS_utils.c
   @date      Sep 17, 2016
   @author    Daniel Siegel, Philipp Moesta
   @desc 
   Calls to EOS_Omni routines required by con2prim.
   @enddesc 
 @@*/
 

#include "con2prim.h"

#if !STANDALONE
  #include "cctk.h"
  #include "cctk_Arguments.h"
  #include "cctk_Parameters.h"
  #include "cctk_Functions.h"
#endif


extern struct c2p_steer c2p;


void EOS_press(int eoskey, int keytemp, double rho, double * eps, 
        double * temp, double ye, double * prs, int * keyerr, int * nEOS_calls){

#if STANDALONE
  const int eos_key = eoskey;
  const int havetemp = keytemp;
  const double prec = c2p.eos_prec;
  int keyerr_loc;
  double xrho, xeps, xtemp, xye, xprs;
#else
  const CCTK_INT eos_key = eoskey;
  const CCTK_INT havetemp = keytemp;
  const CCTK_REAL prec = c2p.eos_prec;
  const CCTK_INT npoints = 1;
  CCTK_INT keyerr_loc;
  CCTK_INT anyerr = 0;
  CCTK_REAL xrho, xeps, xtemp, xye, xprs;
#endif

  int nEOScalls = 0;
  keyerr_loc = 0;

  xrho = rho;
  xtemp = *temp;
  xeps = *eps;
  xye = ye;
  xprs = 0.0;
  
#if STANDALONE
  EOS_Omni_press(eos_key,havetemp,prec,xrho,&xeps,&xtemp,xye,&xprs,&keyerr_loc,
    &nEOScalls);
#else
  EOS_Omni_press(eos_key,havetemp,prec,npoints,&xrho,&xeps,&xtemp,&xye,&xprs,
    &keyerr_loc,&anyerr);
#endif

  *eps  = xeps;
  *temp = xtemp;
  *prs  = xprs;
  *keyerr = keyerr_loc;
  *nEOS_calls = nEOScalls;
}


void EOS_press_ent_abar(int eoskey, int keytemp, double rho, double * eps, 
        double * temp, double ye, double * prs, double * ent, double * abar, 
        int * keyerr, int * nEOS_calls){

#if STANDALONE
  const int eos_key = eoskey;
  const int havetemp = keytemp;
  const double prec = c2p.eos_prec;
  int keyerr_loc;
  double xrho, xeps, xtemp, xye, xprs, xent, xabar;
#else
  const CCTK_INT eos_key = eoskey;
  const CCTK_INT havetemp = keytemp;
  const CCTK_REAL prec = c2p.eos_prec;
  const CCTK_INT npoints = 1;
  CCTK_INT keyerr_loc;
  CCTK_INT anyerr = 0;
  CCTK_REAL xrho, xeps, xtemp, xye, xprs, xent, xabar;
#endif

  int nEOScalls = 0;
  keyerr_loc = 0;

  xrho = rho;
  xtemp = *temp;
  xeps = *eps;
  xye = ye;
  xprs = 0.0;
  xent = 0.0;
  xabar = 0.0;
  
#if STANDALONE
  EOS_Omni_press_ent_abar(eos_key,havetemp,prec,xrho,&xeps,&xtemp,xye,&xprs,&xent,&xabar,
    &keyerr_loc,&nEOScalls);
#else
  EOS_Omni_press_ent_abar(eos_key,havetemp,prec,npoints,&xrho,&xeps,&xtemp,&xye,&xprs,
    &xent,&xabar,&keyerr_loc,&anyerr);
#endif

  *eps  = xeps;
  *temp = xtemp;
  *prs  = xprs;
  *ent  = xent;
  *abar = xabar;
  *keyerr = keyerr_loc;
  *nEOS_calls = nEOScalls;
}


void EOS_eps_from_press(int eoskey, double rho, double * eps, 
        double * temp, double ye, double prs, int * keyerr){

#if STANDALONE
  const int eos_key = eoskey;
  const int havetemp = 0;
  const double prec = c2p.eos_prec;
  int keyerr_loc;
  double xrho, xeps, xtemp, xye, xprs;
#else
  const CCTK_INT eos_key = eoskey;
  const CCTK_INT havetemp = 0;
  const CCTK_REAL prec = c2p.eos_prec;
  const CCTK_INT npoints = 1;
  CCTK_INT keyerr_loc;
  CCTK_INT anyerr = 0;
  CCTK_REAL xrho, xeps, xtemp, xye, xprs;
#endif

  keyerr_loc = 0;

  xrho = rho;
  xtemp = 0.0;
  xeps = 0.0;
  xye = ye;
  xprs = prs;
  
#if STANDALONE
  EOS_Omni_EpsFromPress(eos_key,havetemp,prec,
                   xrho,&xeps,&xtemp,xye,xprs,&keyerr_loc);
#else
  EOS_Omni_EpsFromPress(eos_key,havetemp,prec,npoints,
                   &xrho,&xeps,&xtemp,&xye,&xprs,&xeps,&keyerr_loc,&anyerr);
#endif

  *eps  = xeps;
  *temp = xtemp;
  *keyerr = keyerr_loc;

}


void EOS_EP_dEdr_dEdt_dPdr_dPdt(double * x, const double * con, 
          double * Eprim, double * Pprim, double * dEdrho, double * dEdt, double * dPdrho, 
          double * dPdt, int stepsize){
  
  /*  Compute partial derivatives of specific internal energy and pressure with respect
   *  to density and temperature, based on primitives computed from Newton-Raphson state
   *  vector x and conservatives
   */
  
  const double W = x[0];
  const double Z = x[1];
  const double T = x[2];

#if STANDALONE
  const int eos_key = c2p.eoskey; 
  const int keytemp = 1;
  const int step_size = stepsize;
  const double prec = c2p.eos_prec;
  int keyerr;
  double xrho,xtemp,xeps,xye,xprs,xdedrho,xdpdrho,xdedt,xdpdt;
#else
  const CCTK_INT eos_key = c2p.eoskey;
  const CCTK_INT keytemp = 1;
  const CCTK_INT step_size = stepsize;
  const CCTK_REAL prec = c2p.eos_prec;
  CCTK_INT keyerr;
  CCTK_REAL xrho,xtemp,xeps,xye,xprs,xdedrho,xdpdrho,xdedt,xdpdt;
#endif  
  
  keyerr = 0;
  
  xrho = con[D]/W;
  xtemp = T;
  xeps = 0.0;
  xye = con[YE]/con[D];
  xprs = 0.0;
  
  xdedrho = 0.0;
  xdpdrho = 0.0;
  xdedt = 0.0;
  xdpdt = 0.0;

  
#if STANDALONE  
  EOS_Omni_dpdrho_dpdt_dedrho_dedt(eos_key, keytemp, xrho, &xtemp, xye, &xeps, &xprs, &xdpdrho, 
        &xdpdt, &xdedrho, &xdedt, &keyerr, step_size);
#else
  EOS_Omni_dpdrho_dpdt_dedrho_dedt(eos_key, keytemp, xrho, &xtemp, xye, &xeps, &xprs, &xdpdrho, 
        &xdpdt, &xdedrho, &xdedt, &keyerr, step_size);
#endif

#if DEBUG
  if (keyerr != 0) printf("EOS_EP_dEdr_dEdt_dPdr_dPdt: keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g p: %g\n",keyerr,xtemp,xrho,xeps,xye,xprs);
#endif    
    
  *dEdrho = xdedrho;
  *dEdt   = xdedt;
  *dPdrho = xdpdrho;
  *dPdt   = xdpdt;
  *Eprim  = xeps;
  *Pprim  = xprs;
  
}


void EOS_EP_dEdr_dEdt_dPdr_dPdt_2D(const double rho2D, const double temp2D, const double * con, 
          double * Eprim, double * Pprim, double * dEdrho, double * dEdt, double * dPdrho, 
          double * dPdt, int stepsize){
  
  /*  Compute partial derivatives of specific internal energy and pressure with respect
   *  to density and temperature, based on primitives computed from Newton-Raphson state
   *  vector x and conservatives
   */

#if STANDALONE
  const int eos_key = c2p.eoskey; 
  const int keytemp = 1;
  const int step_size = stepsize;
  const double prec = c2p.eos_prec;
  int keyerr;
  double xrho,xtemp,xeps,xye,xprs,xdedrho,xdpdrho,xdedt,xdpdt;
#else
  const CCTK_INT eos_key = c2p.eoskey;
  const CCTK_INT keytemp = 1;
  const CCTK_INT step_size = stepsize;
  const CCTK_REAL prec = c2p.eos_prec;
  CCTK_INT keyerr;
  CCTK_REAL xrho,xtemp,xeps,xye,xprs,xdedrho,xdpdrho,xdedt,xdpdt;
#endif  
  
  keyerr = 0;
  
  xrho = rho2D; //con[D]/x[0];
  xtemp = temp2D; //x[1];
  xeps = 0.0;
  xye = con[YE]/con[D];
  xprs = 0.0;
  
  xdedrho = 0.0;
  xdpdrho = 0.0;
  xdedt = 0.0;
  xdpdt = 0.0;

  
#if STANDALONE  
  EOS_Omni_dpdrho_dpdt_dedrho_dedt(eos_key, keytemp, xrho, &xtemp, xye, &xeps, &xprs, &xdpdrho, 
        &xdpdt, &xdedrho, &xdedt, &keyerr, step_size);
#else
  EOS_Omni_dpdrho_dpdt_dedrho_dedt(eos_key, keytemp, xrho, &xtemp, xye, &xeps, &xprs, &xdpdrho, 
        &xdpdt, &xdedrho, &xdedt, &keyerr, step_size);
#endif

#if DEBUG
  if (keyerr != 0) printf("EOS_EP_dEdr_dEdt_dPdr_dPdt: keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g p: %g\n",keyerr,xtemp,xrho,xeps,xye,xprs);
#endif    

  *dEdrho = xdedrho;
  *dEdt   = xdedt;
  *dPdrho = xdpdrho;
  *dPdt   = xdpdt;
  *Eprim  = xeps;
  *Pprim  = xprs;
  
}


void EOS_P_from_hrho_dPdrho_dPdeps(const double rho, const double enth, const double * con, 
          double * temp_guess, double * eps, double * press, double * dPdrho, 
          double * dPdeps, double * entr, double * abar, int * nEOS_calls){
  
  /*  Perform inversion from specific enthalpy to temperature and compute partial derivatives 
   *  of pressure wrt density and specific internal energy
   */
  
#if STANDALONE
  const int eos_key = c2p.eoskey;
  const int step_size = 1;
  const double prec = c2p.eos_prec;
  int keyerr, xnEOS_calls;
  double xtemp,xeps,xprs,xdpdrho,xdpdeps,xentr,xabar;
  const double xenth = enth;
  const double xrho = rho;
  const double xye = con[YE]/con[D];
#else
  const CCTK_INT eos_key = c2p.eoskey;
  const CCTK_INT step_size = 1;
  const CCTK_REAL prec = c2p.eos_prec;
  const CCTK_REAL xrho = rho;
  const CCTK_REAL xenth = enth;
  const CCTK_REAL xye = con[YE]/con[D];
  CCTK_INT keyerr;
  CCTK_REAL xtemp,xeps,xprs,xdpdrho,xdpdeps,xentr,xabar;
#endif  
  
  keyerr = 0;
  xtemp = *temp_guess;
  xeps = 0.0;
  xprs = 0.0;
  xdpdeps = 0.0;
  xdpdrho = 0.0;
  xentr = 0.0;
  xabar = 0.0;
  
#if STANDALONE  
  EOS_Omni_dpdrho_dpdeps_hinv(eos_key, xrho, xenth, &xtemp, xye, &xeps, &xprs, &xdpdrho, 
        &xdpdeps, &xentr, &xabar, &keyerr, step_size, &xnEOS_calls);
  *nEOS_calls = xnEOS_calls;      
#else
  printf("%s\n", "EOS_press_from_rhoenthalpy not yet implemented for evolution code!");
  exit(1);
#endif

#if DEBUG
  if (keyerr != 0) printf("EOS_P_from_hrho_dPdr_dPdeps: keyerr: %d, Temp_guess: %g Rho: %g enthalpy: %g Ye: %g p: %g\n",keyerr,xtemp,xrho,xenth,xye,xprs);
#endif    
    
  *temp_guess = xtemp;
  *dPdeps = xdpdeps;
  *dPdrho = xdpdrho;
  *press  = xprs;
  *eps = xeps;
  *entr = xentr;
  *abar = xabar;
  
}


void EOS_press_from_rhoenthalpy(int eoskey, int keytemp, double prec, double rho, double * eps, 
    double * temp, double ye, double * press, double * enth, int * anyerr, int * keyerr, int * nEOScalls)
{  
  /*  Perform inversion from specific enthalpy to temperature and compute pressure
   *  and specific internal energy
   */
  
#if STANDALONE
  const int eos_key = c2p.eoskey;
  const int key_temp = keytemp;
  const double preci = prec;
  int xkeyerr, xanyerr, xnEOScalls;
  double xtemp,xeps,xprs;
  double xenth = *enth;
  const double xrho = rho;
  const double xye = ye;
#else
  const CCTK_INT eos_key = c2p.eoskey;
  const CCTK_INT key_temp = keytemp;
  const CCTK_REAL preci = prec;
  const CCTK_REAL xrho = rho;
  const CCTK_REAL xye = ye;
  CCTK_REAL xenth = *enth;
  CCTK_INT xkeyerr, xanyerr;
  CCTK_REAL xtemp,xeps,xprs;
#endif  
  
  xkeyerr = 0;
  xanyerr = 0;
  xtemp = *temp;
  xeps = 0.0;
  xprs = 0.0;
  
#if STANDALONE
  EOS_Omni_press_from_rhoenthalpy(eos_key, key_temp, preci, xrho, &xeps, 
    &xtemp, xye, &xprs, &xenth, &xanyerr, &xkeyerr, &xnEOScalls);
  *nEOScalls = xnEOScalls;      
#else
  printf("%s\n", "EOS_press_from_rhoenthalpy not yet implemented for evolution code!");
  exit(1);
#endif

#if DEBUG
  if (keyerr != 0) printf("EOS_press_from_rhoenthalpy: keyerr: %d, Temp_guess: %g Rho: %g enthalpy: %g Ye: %g p: %g\n",keyerr,xtemp,xrho,xenth,xye,xprs);
#endif    
    
  *temp = xtemp;
  *press  = xprs;
  *eps = xeps;
  *keyerr = xkeyerr;
  *anyerr = xanyerr;

}