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

extern double EOS_const[5];
extern void helm_eos_m_kt1_dpdrho_dpdt_dedrho_dedt_point_(double * xrho, 
    double * xtemp, double * xye, double * xeps, double * xpress, double * xdpdrho,
    double * xdpdt, double * xdedrho, double * xdedt, int * xanyerr);
extern void helm_eos_kt1_press_eps_point_(double * xrho, double * xtemp, 
    double * xye, double * xeps, double * xprs, int * anyerr); 
extern void helm_eos_kt0_press_temp_point_(double * xrho, double * xtemp, 
    double * xye, double * xeps, double * xprs, int * anyerr, int * nEOScalls);
extern void helm_eos_kt1_press_eps_ent_abar_point_(double * xrho, double * xtemp, 
    double * xye, double * xeps, double * xprs, double * ent, double * abar, int * anyerr); 
extern void helm_eos_kt0_press_temp_ent_abar_point_(double * xrho, double * xtemp, 
    double * xye, double * xeps, double * xprs, double * ent, double * abar, int * anyerr,
    int * nEOScalls);
extern void helm_eos_m_kt0_dpressdrho_dpressdeps_(int npoints, double * xrho, double * xtemp, 
    double * xye, double * xeps, double * xdpdrho, double * xdpdeps, int * inv_err, int * anyerr, 
    int * xnEOScalls);
extern void helm_eos_m_kt0_dpressdrho_dpressdeps_hinv_(int * npoints, double * xrho, 
    double * xtemp, double * xenth, double * xye, double * xeps, double * xprs, 
    double * xdpdrho, double * xdpdeps, double * xentr, double * xabar, int * inv_err, int * anyerr, int * xnEOScalls);
extern void helm_eos_m_kt1_dpressdrho_dpressdeps_dpressdt_depsdt_(int * npoints,double * xrho,
    double * xtemp, double * xye, double * xeps, double * xprs, double * xdpdrho, double * xdpdeps,
    double * xdpdt, double * xdepsdt, int * keyerr, int * anyerr);

void nuc_eos_P_from_Enthalpy(double rho, double *temp, double ye,
		     double *enr, double *enr2, double *press,
		     int keytemp,
		     int *keyerr,double rfeps,
         int * xnEOScalls); 

void EOS_Omni_press(int eoskey, int keytemp, double prec, double rho, double * eps, 
    double * temp, double ye, double * press, int * keyerr, int * nEOScalls){

  *nEOScalls = 1; //otherwise reset below

  if(eoskey==1){ // POLYTROPIC
    
    double gammapoly = EOS_const[0];
    double Kpoly = EOS_const[1];
    
    *press = Kpoly*pow(rho,gammapoly);
    
    if(keytemp==1){
      *eps = *press / ((gammapoly - 1.0) * rho);
    }
  }
  else if(eoskey==2){ // IDEAL GAS
    const double gammaideal = EOS_const[0];
    const double eps_over_T_IG = (kBerg / (amu * (gammaideal - 1.0))) * eps_gf * inv_temp_gf;
    
    if(keytemp==1){
      *eps = *temp * eps_over_T_IG;
    }
    else if (keytemp==0){
      *temp  = *eps / eps_over_T_IG;
    }
    else{
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    }

    *press = *eps * rho * (gammaideal - 1.0); 

  }
  else if(eoskey==3){ // Hybrid EOS

    double gamma1 = EOS_const[0];
    double K1 = EOS_const[1];
    double gamma2 = EOS_const[2];
    double gammath = EOS_const[3];
    double rho_nuc = EOS_const[4];
    
    // Equations 8-9 from Janka et al. 1993
    double E1 = K1/(gamma1-1.0);
    double E2 = (gamma1-1.0)/(gamma2-1.0)*E1*pow(rho_nuc,(gamma1-gamma2)); 
    double K2 = (gamma2-1.0)*E2;
    double E3 = (gamma2-gamma1)/(gamma2-1.0)*E1*pow(rho_nuc,(gamma1-1.0));
    
    double Kx, Ex, Gx, Ex3;
    double eps_p, eps_th,p_th,p_poly;
    
    if(rho<rho_nuc){
      Kx = K1;
      Ex = E1;
      Gx = gamma1;
      Ex3 = 0.0;
    }
    else{
      Kx = K2;
      Ex = E2;
      Gx = gamma2;
      Ex3 = E3;
    }
    
    // Equation (7) from Janka et al. 1993
    eps_p = Ex*pow(rho,Gx) + Ex3*rho;
    
    eps_th = *eps * rho - eps_p;
    
    // Equation (11) from Janka et al. 1993
    p_th = (gammath-1.0)*eps_th;
    
    p_th = p_th > 0.0 ? p_th : 0.0;
   
    p_poly = Kx*pow(rho,Gx);
    
    *press = p_poly + p_th;
  }
  else if(eoskey==4){ //nuc eos

    int xnEOScalls = 0;
    int npoints = 1;
    int anyerr = 0;

    //printf("EOS_Omni Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    if(keytemp==1) {
      nuc_eos_m_kt1_press_eps(&npoints,&rho,temp,&ye,eps,press,keyerr,&anyerr);
    } else if (keytemp==0) {
      nuc_eos_m_kt0_press(&npoints,&rho,temp,&ye,eps,press,&prec,keyerr,&anyerr,&xnEOScalls);
      *nEOScalls = xnEOScalls;
    } else {
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    }

    if (*keyerr != 0) printf("EOS_Omni Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
  
  }
  else if(eoskey==8){ //Helmholtz EOS
    
    double xrho = rho;
    double xye = ye;
    double xtemp = *temp;
    double xeps = *eps;
    double xprs = 0.0;
    int anyerr = 0;

    if(keytemp==1){
      helm_eos_kt1_press_eps_point_(&xrho, &xtemp, &xye, &xeps, &xprs, &anyerr);
      *eps = xeps;
    } else if (keytemp==0) {
      int xnEOScalls = 0;
      helm_eos_kt0_press_temp_point_(&xrho, &xtemp, &xye, &xeps, &xprs, &anyerr, &xnEOScalls);
      *temp = xtemp;
      *nEOScalls = xnEOScalls;
    } else {
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    }

    if (anyerr != 0) printf("EOS_Omni ERROR: Temp: %g Rho: %g Eps: %g Ye: %g\n",*temp,rho,*eps,ye);
    
    *press = xprs;
  
  }
  else {
    printf("This EOS identifier is not supported! Stopping the code...\n");
    exit(EXIT_FAILURE);
  }
}



void EOS_Omni_press_ent_abar(int eoskey, int keytemp, double prec, double rho, double * eps, 
    double * temp, double ye, double * press, double * ent, double * abar, int * keyerr, 
    int * nEOScalls){

  *nEOScalls = 1; //otherwise reset below
  
  if(eoskey==4){ //nuc eos

    int xnEOScalls = 0;
    int npoints = 1;
    int anyerr = 0;

    //printf("EOS_Omni Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    if(keytemp==1) {
      nuc_eos_m_kt1_press_eps(&npoints,&rho,temp,&ye,eps,press,keyerr,&anyerr);
    } else if (keytemp==0) {
      nuc_eos_m_kt0_press(&npoints,&rho,temp,&ye,eps,press,&prec,keyerr,&anyerr,&xnEOScalls);
      *nEOScalls = xnEOScalls;
    } else {
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    }

    if (*keyerr != 0) printf("EOS_Omni Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    
    // Note: entropy and abar are not yet returned
    *ent = 0.0; // set to zero for now
    *abar = 0.0; // set to zero for now
  }
  else if(eoskey==8){ //Helmholtz EOS
    
    double xrho = rho;
    double xye = ye;
    double xtemp = *temp;
    double xeps = *eps;
    double xprs = 0.0;
    double xent = 0.0;
    double xabar = 0.0;
    int anyerr = 0;

    if(keytemp==1){
      helm_eos_kt1_press_eps_ent_abar_point_(&xrho, &xtemp, &xye, &xeps, &xprs, &xent, &xabar, &anyerr);
      *eps = xeps;
    } else if (keytemp==0) {
      int xnEOScalls = 1;
      helm_eos_kt0_press_temp_ent_abar_point_(&xrho, &xtemp, &xye, &xeps, &xprs, &xent, &xabar, &anyerr, &xnEOScalls);
      *temp = xtemp;
      *nEOScalls = xnEOScalls;
    } else {
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    } 

    if (anyerr != 0) printf("EOS_Omni ERROR: Temp: %g Rho: %g Eps: %g Ye: %g\n",*temp,rho,*eps,ye);
    
    *press = xprs;
    *ent = xent;
    *abar = xabar;
  
  }
  else {
    printf("This EOS identifier is not supported! Stopping the code...\n");
    exit(EXIT_FAILURE);
  }
}


void EOS_Omni_EpsFromPress(int eoskey, int keytemp, double prec, double rho, double * eps, 
    double * temp, double ye, double press, int * keyerr)
{
    
  if(eoskey==1){ // POLYTROPIC
    
    const double gammapoly = EOS_const[0];
    const double Kpoly = EOS_const[1];
    *eps = press / ((gammapoly - 1.0) * rho);
  }
  else if(eoskey==2){ // IDEAL GAS
    const double gammaideal = EOS_const[0];
    const double eps_over_T_IG = (kBerg / (amu * (gammaideal - 1.0))) * eps_gf * inv_temp_gf;
    
    *eps = press / ((gammaideal - 1.0) * rho);
    
    if (keytemp==0){
      *temp  = *eps / eps_over_T_IG;
    }
  } else {
    printf("This EOS identifier is not supported! Stopping the code...\n");
    exit(EXIT_FAILURE);
  }
}

void EOS_Omni_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){

  if(eoskey==1){ // POLYTROPIC
    
    double gammapoly = EOS_const[0];
    double Kpoly = EOS_const[1];
//    printf("gammapoly=%e, Kpoly=%e, rho=%e", gammapoly, Kpoly, rho); 
    *press = Kpoly*pow(rho,gammapoly);
    *nEOScalls = 1;
    }
  else if(eoskey==2){ //IDEAL GAS
    double gammaideal = EOS_const[0];
    *press = (*enth - rho) * (gammaideal-1.0)/gammaideal;
    *eps = *press/(rho * (gammaideal-1.0));
    *nEOScalls = 1;
  }
  else if(eoskey==4){ //nuc eos

    double enthm1;
    const double prec = 1.0e-10;
    
    enthm1 = (*enth/rho-1.0);
    
#if DEBUG
    printf("EOS_Omni before call; Rho: %g Eps: %g Ye: %g, enth = %g, enthm1 = %g\n",rho,*eps,ye,*enth, enthm1);
#endif    

    *nEOScalls = 0;
    int xnEOScalls;
    nuc_eos_P_from_Enthalpy(rho, temp, ye, &enthm1, eps, press,
		    keytemp, keyerr, prec, &xnEOScalls);
#if DEBUG
    printf("update EOSOmni: nEOScalls = %i\n", xnEOScalls);
#endif    

    *nEOScalls += xnEOScalls;

    if (*keyerr != 0) printf("EOS_Omni Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    
  }
  else if (eoskey==8){ //Helmholtz EOS
    double xrho = rho;
    double xenth = (*enth/rho-1.0); // Helmholtz EOS takes non-relativistic enthalpy
    double xtemp = *temp;
    double xye = ye;
    double xeps = 0.0;
    double xprs = 0.0;

    double xdpdrho = 0.0;
    double xdpdeps = 0.0;
    double xentr = 0.0;
    double xabar = 0.0;
    
    int anyerr = 0;
    int inv_err = 0;
    int npoints = 1;
    int xnEOScalls = 0;
    
    helm_eos_m_kt0_dpressdrho_dpressdeps_hinv_(&npoints, &xrho, &xtemp, &xenth, &xye, &xeps, &xprs,
        &xdpdrho, &xdpdeps, &xentr, &xabar, &inv_err, &anyerr, &xnEOScalls);
    
    if (anyerr != 0) printf("EOS_Omni_dEp_hinv ERROR: Temp_guess: %g Rho: %g Eps: %g Ye: %g\n",*temp,rho,*enth,ye);
    
    *nEOScalls = xnEOScalls;
    *eps = xeps;
    *temp = xtemp;
    *press = xprs;
  }
  else{
    printf("This EOS identifier is not supported! Stopping the code...\n");
    exit(EXIT_FAILURE);
  }

}
void EOS_Omni_eps_from_rhopress(int eoskey, int keytemp, double prec, int n, double rho, double * eps, 
    double * temp, double ye, double * press,int * anyerr, int * keyerr){
    if(eoskey==1){
    if(keytemp==1){
      double gammapoly = EOS_const[0];
      double Kpoly = EOS_const[1];
      *eps = *press / ((gammapoly - 1.0) * rho);
    }
    }
    if(eoskey==2){
      double gammaideal = EOS_const[0];
      *eps = *press / ((gammaideal - 1.0) * rho);
    }
  }

void EOS_Omni_enthalpy_from_rhopress(int eoskey, int keytemp, double prec, int n, double rho, double * eps, 
    double * temp, double ye, double * press, double * ent, int * anyerr, int * keyerr){
    if(eoskey==1){
    if(keytemp==1){
      double gammapoly = EOS_const[0];
      double Kpoly = EOS_const[1];
      *eps = *press / ((gammapoly - 1.0) * rho);
      *ent = 1.0 + *eps + *press/rho;
    }
    }
  }


void EOS_Omni_get_eos_limits(int eoskey, double * rho_min, double * rho_max, double * t_min, 
        double * t_max, double * ye_min, double * ye_max, double * eps_min, double * eps_max, double * press_min, double * press_max)
{ 
  // default settings
  *t_min = 0.0;
  *t_max = 200.0;
  *ye_min = 0.0;
  *ye_max = 1.0;
  *rho_min = 0.0;
  *rho_max = 1.0e16;
  *eps_min = 0.0;
  *eps_max = 1.0e50;    
      
  if(eoskey==4){ // nuc EOS
    *rho_min = eos_rhomin;
    *rho_max = eos_rhomax;
    *t_min   = eos_tempmin;
    *t_max   = eos_tempmax;
    *ye_min  = eos_yemin;
    *ye_max  = eos_yemax;
    *eps_min = eos_epsmin;
    *eps_max = eos_epsmax;
    *press_min = eos_pressmin;
    *press_max = eos_pressmax;

#if DEBUG    
    printf("rho min, max: %g \t %g\n", *rho_min, *rho_max);
    printf("temp min, max: %g \t %g\n", *t_min, *t_max);
    printf("ye min, max: %g \t %g\n", *ye_min, *ye_max);
    printf("eps min, max: %g \t %g\n", *eps_min, *eps_max);
    printf("press min, max: %g \t %g\n", *press_min, *press_max);
#endif

  } else if (eoskey==8){ //Helmholtz EOS      

    *t_min = 8.62e-7; // 8.62e-7 = 1e4K: lower limit of HelmEOS table
    *eps_min = -0.0075941539; //HelmEOS: -0.007594153911 = - Qalpha/(4*m_u)
  }
}

void EOS_Omni_dpdrho_dpdt_dedrho_dedt(int eoskey, int keytemp, double rho, double * temp, 
       double ye, double * eps, double * press, double * dpdrho, double * dpdt, double * dedrho, 
       double * dedt, int * keyerr, int stepsize)
{

  if(eoskey==1){ // POLYTROPIC
        
    const double gammapoly = EOS_const[0];
    const double Kpoly = EOS_const[1];
    
    *press = Kpoly*pow(rho,gammapoly);
    
    if(keytemp==1){
      *eps = *press / ((gammapoly - 1.0) * rho);      
    }
        
    *dpdrho = Kpoly * gammapoly * pow(rho,gammapoly-1.0);
    *dpdt = 0.0;
    *dedrho = Kpoly * pow(rho,gammapoly - 2.0);
    *dedt = 0.0;
  }
  else if(eoskey==2){ //IDEAL GAS
    
    const double gammaideal = EOS_const[0];
    const double eps_over_T_IG = (kBerg / (amu * (gammaideal - 1.0))) * eps_gf * inv_temp_gf;
    
    if(keytemp==1){
      *eps = *temp * eps_over_T_IG;
    }
    else if (keytemp==0){
      *temp  = *eps / eps_over_T_IG;
    }
    else{
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    }
            
    *press = *eps * rho * (gammaideal - 1.0);
    *dpdrho = *eps * (gammaideal-1.0);
    *dpdt   = rho * (kBerg / amu) * eps_gf * inv_temp_gf;
    *dedrho = 0.0;
    *dedt   = eps_over_T_IG;
        
  }      
  else if(eoskey==3){ // Hybrid EOS      
        
    printf("This EOS is not yet fully implemented! Stopping the code...\n");
    exit(EXIT_FAILURE);      
  
  }
  else if(eoskey==4){ // nuc EOS
  
    double xtemp,xeps,xprs;
    const double prec = 1.0e-10;
    double xeps2 = 0.0; 
    int anyerr = 0; 

    xtemp = *temp;
    xeps = *eps * inv_eps_gf;
    xprs = 0.0;

    double xdpdrho = 0.0;
    double xdedt = 0.0;
    double xdpdt = 0.0;
    double xdedrho = 0.0;
    int npoints =1;  

    nuc_eos_m_kt1_dpdrhot_dedrhot(&npoints,&rho,temp,&ye,eps,&xeps2,press,
                    &xdpdrho,&xdpdt,&xdedrho,&xdedt,&prec,keyerr,&anyerr);

    if (*keyerr != 0) printf("EOS_Omni_dEp Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    
    *dedrho = (xeps2 / rho) * xdedrho;
    *dedt   = (xeps2 / *temp) * xdedt; 
    *dpdrho = (*press / rho) * xdpdrho;
    *dpdt   = (*press / *temp) * xdpdt; 

  }

  else if(eoskey==8){ // Helmholtz EOS      
    
    const double prec = 1.0e-10;
    double xrho = rho;
    double xtemp = *temp;
    double xye = ye;
    double xeps = *eps;
    double xprs = 0.0;

    double xdedrho = 0.0;
    double xdpdrho = 0.0;
    double xdedt = 0.0;
    double xdpdt = 0.0;
    int anyerr = 0;
    
    if(keytemp==1){  
      helm_eos_m_kt1_dpdrho_dpdt_dedrho_dedt_point_(&xrho, &xtemp, &xye, &xeps, &xprs, 
               &xdpdrho,&xdpdt, &xdedrho, &xdedt, &anyerr);
      *eps = xeps;
    }
    else{
      printf("This keytemp is not suppported yet! Stopping the code...");
      exit(EXIT_FAILURE);
    }
    
    if (anyerr != 0) printf("EOS_Omni_dEp ERROR: Temp: %g Rho: %g Eps: %g Ye: %g\n",*temp,rho,*eps,ye);
    
    *press = xprs;
    *dedrho = xdedrho;
    *dedt   = xdedt;
    *dpdrho = xdpdrho;
    *dpdt   = xdpdt;
       
  } else {     
    printf("eoskey %d is not suppported yet! Stopping the code...\n", eoskey);
    exit(EXIT_FAILURE);
  }

}

void EOS_Omni_dpdrho_dpdeps(const int eoskey, const int keytemp, const double rho, double * temp, 
       const double ye, double * eps, double * press, double * dpdrho, double * dpdeps, int * keyerr,
       int stepsize, int * nEOScalls) {

  if(eoskey==1){ // POLYTROPIC
        
    const double gammapoly = EOS_const[0];
    const double Kpoly = EOS_const[1];
    
    *press = Kpoly*pow(rho,gammapoly);
    
    if(keytemp==1){
      *eps = *press / ((gammapoly - 1.0) * rho);      
    }
    
    *nEOScalls = 1;   
    *dpdrho = Kpoly * gammapoly * pow(rho,gammapoly-1.0);
    *dpdeps = 0.0;
  }
  else if(eoskey==2){ //IDEAL GAS
    
    const double gammaideal = EOS_const[0];
    const double eps_over_T_IG = (kBerg / (amu * (gammaideal - 1.0))) * eps_gf * inv_temp_gf;
    
    if(keytemp==1){
      *eps = *temp * eps_over_T_IG;
    }
    else if (keytemp==0){
      *temp  = *eps / eps_over_T_IG;
    }
    else{
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    }
    
    *nEOScalls = 1;        
    *press = *eps * rho * (gammaideal - 1.0);
    *dpdrho = *eps * (gammaideal-1.0);
    *dpdeps = rho * (gammaideal-1.0);
  }      
  else if(eoskey==3){ // Hybrid EOS      
        
    printf("This EOS is not yet fully implemented! Stopping the code...\n");
    exit(EXIT_FAILURE);      
  
  }
  else if(eoskey==4){ // nuc EOS
  
    printf("This EOS is not yet fully implemented! Stopping the code...\n");
    exit(EXIT_FAILURE);      

  }
  else if(eoskey==8){ // Helmholtz EOS      
    
    //printf("This EOS is not yet fully implemented! Stopping the code...\n");
    //exit(EXIT_FAILURE);
    const double prec = 1.0e-10;
    double xrho = rho;
    double xtemp = *temp;
    double xye = ye;
    double xeps = *eps;
    double xprs = 0.0;

    double xdpdrho = 0.0;
    double xdpdeps = 0.0;
    
    int anyerr = 0;
    int inv_err = 0;
    int npoints = 1;
    int xnEOScalls = 0;
    
    if(keytemp==1){  
      printf("This keytemp is not suppported yet! Stopping the code...");
      exit(EXIT_FAILURE);
    }
    else if (keytemp==1){
      helm_eos_m_kt0_dpressdrho_dpressdeps_(npoints, &xrho, &xtemp, &xye, &xeps, &xdpdrho, &xdpdeps,
        &inv_err, &anyerr, &xnEOScalls);
      *eps = xeps;  
    } else {
      printf("This keytemp is not suppported yet! Stopping the code...");
      exit(EXIT_FAILURE);
    }
    
    if (anyerr != 0) printf("EOS_Omni_dEp ERROR: Temp: %g Rho: %g Eps: %g Ye: %g\n",*temp,rho,*eps,ye);
    
    *nEOScalls = xnEOScalls;
    *press = xprs;
    *dpdrho = xdpdrho;
    *dpdeps = xdpdeps;

  } else {     
    printf("eoskey %d is not suppported yet! Stopping the code...\n", eoskey);
    exit(EXIT_FAILURE);
  }

 }


void EOS_Omni_dpdrho_dpdeps_hinv(const int eoskey, const double rho, const double enth, double * temp, 
       const double ye, double * eps, double * press, double * dpdrho, double * dpdeps, double * entr, 
       double * abar, int * keyerr, int stepsize, int * nEOScalls) {

  if(eoskey==1){ // POLYTROPIC
        
    const double gammapoly = EOS_const[0];
    const double Kpoly = EOS_const[1];
    
    *nEOScalls = 1;  
    *press = Kpoly*pow(rho,gammapoly);
    *eps = *press / ((gammapoly - 1.0) * rho); 
    *dpdrho = Kpoly * gammapoly * pow(rho,gammapoly-1.0);
    *dpdeps = 0.0;
  }
  else if(eoskey==2){ //IDEAL GAS
    
    const double gammaideal = EOS_const[0];
    const double eps_over_T_IG = (kBerg / (amu * (gammaideal - 1.0))) * eps_gf * inv_temp_gf;
    
    *nEOScalls = 1;
    *eps = (enth - 1.0)/gammaideal;
    *temp  = *eps / eps_over_T_IG; 
    *press = *eps * rho * (gammaideal - 1.0);
    *dpdrho = *eps * (gammaideal-1.0);
    *dpdeps = rho * (gammaideal-1.0);
  }      
  else if(eoskey==3){ // Hybrid EOS      
        
    printf("This EOS is not yet fully implemented! Stopping the code...\n");
    exit(EXIT_FAILURE);      
  
  }
  else if(eoskey==4){ // nuc EOS

    double entm1;
    const double prec = 1.0e-10;
    const int keytemp = 0; 
    entm1 = enth-1.0;
    int xnEOScalls=0;
    *nEOScalls = 0;

    nuc_eos_P_from_Enthalpy(rho, temp, ye, &entm1, eps, press,
		    keytemp, keyerr, prec, &xnEOScalls);
    *nEOScalls += xnEOScalls;

    int anyerr = 0; 
    int npoints =1;  

    nuc_eos_m_kt1_dpdrhoe_dpderho(&npoints,&rho,temp,&ye,eps,
                    dpdrho,dpdeps,&prec,keyerr,&anyerr);
    *nEOScalls += 1;

    if (*keyerr != 0) printf("EOS_Omni_dEp Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    
    // modify EOS call to also output entropy, Abar, set to zero for now
    *entr = 0.0;
    *abar = 0.0;

  }
  else if(eoskey==8){ // Helmholtz EOS      
    
    //const double prec = 1.0e-10;
    double xrho = rho;
    double xenth = enth-1.0; // Helmholtz EOS takes non-relativistic enthalpy
    double xtemp = *temp;
    double xye = ye;
    double xeps = 0.0;
    double xprs = 0.0;
    double xdpdrho = 0.0;
    double xdpdeps = 0.0;
    double xentr = 0.0;
    double xabar = 0.0;

    int anyerr = 0;
    int inv_err = 0;
    int npoints = 1;
    int xnEOScalls = 0;
    
    helm_eos_m_kt0_dpressdrho_dpressdeps_hinv_(&npoints, &xrho, &xtemp, &xenth, &xye, &xeps, &xprs,
        &xdpdrho, &xdpdeps, &xentr, &xabar, &inv_err, &anyerr, &xnEOScalls);
    
    if (anyerr != 0) printf("EOS_Omni_dEp_hinv ERROR: Temp_guess: %g Rho: %g Eps: %g Ye: %g\n",*temp,rho,enth,ye);
    
    *nEOScalls = xnEOScalls;
    *eps = xeps;
    *temp = xtemp;
    *press = xprs;
    *dpdrho = xdpdrho;
    *dpdeps = xdpdeps;
    *entr = xentr;
    *abar = xabar;

  } else {     
    printf("eoskey %d is not suppported yet! Stopping the code...\n", eoskey);
    exit(EXIT_FAILURE);
  }

}


void EOS_Omni_dpdrho_dpdeps_dpdT_depsdT(const int eoskey, const int keytemp, const double rho, double * temp, 
       const double ye, double * eps, double * press, double * dpdrho, double * dpdeps, double * dpdt, double * depsdt, int * keyerr,
       int stepsize, int * nEOScalls) {

  if(eoskey==1){ // POLYTROPIC
        
    double gammapoly = EOS_const[0];
    double Kpoly = EOS_const[1];
    
    *press = Kpoly*pow(rho,gammapoly);
    
    if(keytemp==1){
      *eps = *press / ((gammapoly - 1.0) * rho);      
    }
    
    *nEOScalls = 1;   
    *dpdrho = Kpoly * gammapoly * pow(rho,gammapoly-1.0);
    *dpdeps = 0.0;
    *dpdt = 0.0;
    *depsdt = 0.0;
  }
  else if(eoskey==2){ //IDEAL GAS
    
    const double gammaideal = EOS_const[0];
    const double eps_over_T_IG = (kBerg / (amu * (gammaideal - 1.0))) * eps_gf * inv_temp_gf;
    
    if(keytemp==1){
      *eps = *temp * eps_over_T_IG;
    }
    else if (keytemp==0){
      *temp  = *eps / eps_over_T_IG;
    }
    else{
      printf("This keytemp is not suppported! Stopping the code...\n");
      exit(EXIT_FAILURE);
    }
    
    *nEOScalls = 1;        
    *press = *eps * rho * (gammaideal - 1.0);
    #if DEBUG
    printf("press \t eps \t rho \n");
    printf("%e \t %e \t %e \n", *press, *eps, rho);
    #endif

    *dpdrho = *eps * (gammaideal-1.0);
    *dpdeps = rho * (gammaideal-1.0);
    *dpdt = rho * kBerg/amu * eps_gf * inv_temp_gf;
    *depsdt = eps_over_T_IG;
  }      
  else if(eoskey==3){ // Hybrid EOS      
        
    printf("This EOS is not yet fully implemented! Stopping the code...\n");
    exit(EXIT_FAILURE);      
  
  }
  else if(eoskey==4){ // nuc EOS
  
    const double prec = 1.0e-10;
    int anyerr = 0; 

    double xdpdrho = 0.0;
    double xdpdeps = 0.0;
    double xdedt = 0.0;
    double xdpdt = 0.0;
    double xdedrho = 0.0;
    double xeps2 = 0.0;
    int npoints =1;  

    nuc_eos_m_kt1_dpdrhot_dedrhot_dpdeps(&npoints,&rho,temp,&ye,eps,&xeps2,press,
                    &xdpdrho,&xdpdt,&xdedrho,&xdedt,&xdpdeps,&prec,keyerr,&anyerr);

    if (*keyerr != 0) printf("EOS_Omni_dEp Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    
    *depsdt   = (xeps2 / *temp) * xdedt; 
    *dpdrho = (*press / rho) * xdpdrho;
    *dpdt   = (*press / *temp) * xdpdt; 
    *dpdeps = xdpdeps;
        
    if (*keyerr != 0) printf("EOS_Omni_dEp Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,*eps,ye);
    
  }
  else if(eoskey==8){ // Helmholtz EOS      
    
    const double prec = 1.0e-10;
    double xrho = rho;
    double xtemp = *temp;
    double xye = ye;
    double xeps = *eps;
    double xprs = 0.0;

    double xdpdrho = 0.0;
    double xdpdeps = 0.0;
    double xdpdt   = 0.0;
    double xdepsdt = 0.0;
    
    int anyerr = 0;
    int keyerr = 0;
    int npoints = 1;
    int xnEOScalls = 0;
    
    if(keytemp==1){
      helm_eos_m_kt1_dpressdrho_dpressdeps_dpressdt_depsdt_(&npoints, &xrho, &xtemp, &xye, &xeps, &xprs, &xdpdrho, &xdpdeps,
                                                           &xdpdt, &xdepsdt, &keyerr, &anyerr);
      *eps = xeps;
      xnEOScalls = 1;

    }
    else{
      printf("This keytemp is not suppported yet! Stopping the code...");
      exit(EXIT_FAILURE);
    }
    
    if (anyerr != 0) {
      printf("EOS_Omni_dEp ERROR: Temp: %g Rho: %g Eps: %g Ye: %g\n",*temp,rho,*eps,ye);
    };
    
    *nEOScalls = xnEOScalls;
    *press = xprs;
    *dpdrho = xdpdrho;
    *dpdeps = xdpdeps;
    *dpdt   = xdpdt;
    *depsdt = xdepsdt;

  } else {     
    printf("eoskey %d is not suppported yet! Stopping the code...\n", eoskey);
    exit(EXIT_FAILURE);
  }

}