#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "helpers.hh"

void nuc_eos_m_kt1_dpdrhoe_dpderho(const int *restrict n_in,
				     const double *restrict rho, 
				     double *restrict temp,
				     const double *restrict ye,
				     const double *restrict eps,
				     double *restrict dpdrhoe,
				     double *restrict dpderho,
				     const double *restrict prec,
				     int *restrict keyerr,
				     int *restrict anyerr)
{
  const int n = *n_in;
  int keyerr2;

  *anyerr = 0;

  for(int i=0;i<n;i++) {

    // check if we are fine
    // Note that this code now requires that the
    // temperature guess be within the table bounds
    keyerr2 = checkbounds_kt0_noTcheck(rho[i], ye[i]);
    if(keyerr2 != 0) {
         *anyerr = 1;
         printf("Inside kt1: keyerr: %d ,%g, %g, %g",keyerr2, rho[i], temp[i], ye[i]);
    }
  }

  // Abort if there is any error in checkbounds.
  // This should never happen and the program should abort with
  // a fatal error anyway. No point in doing any further EOS calculations.
  if(*anyerr) return;

  for(int i=0;i<n;i++) {
    const double lr = log(rho[i]);
    const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));

    int idx[8];
    double delx,dely,delz;
    get_interp_spots(lr,lt,ye[i],&delx,&dely,&delz,idx);
     {
	const int iv = 6;
	nuc_eos_C_linterp_one(idx,delx,dely,delz,&(dpdrhoe[i]),iv);
     }
     {
	const int iv = 7;
	nuc_eos_C_linterp_one(idx,delx,dely,delz,&(dpderho[i]),iv);
     }
  }

  return;
}

void nuc_eos_m_kt1_dpdrhot_dedrhot_dpdeps(const int *restrict n_in,
				     const double *restrict rho, 
				     double *restrict temp,
				     const double *restrict ye,
				     double *restrict eps,
				     double *restrict eps2,
				     double *restrict press,
				     double *restrict dpdrho,
				     double *restrict dpdt,
				     double *restrict dedrho,
				     double *restrict dedt,
				     double *restrict dpdeps,
				     const double *restrict prec,
				     int *restrict keyerr,
				     int *restrict anyerr)
{

  const int n = *n_in;
  int keyerr2;

  for(int i=0;i<n;i++) {
    // check if we are fine
    keyerr2 = checkbounds(rho[i], temp[i], ye[i]);
    if(keyerr2 != 0) {
      *anyerr = 1;
      printf("Inside kt1: keyerr: %d ,%g, %g, %g",keyerr2, rho[i], temp[i], ye[i]);
    }
  }

  // Abort if there is any error in checkbounds.
  // This should never happen and the program should abort with
  // a fatal error anyway. No point in doing any further EOS calculations.
  if(*anyerr) return;

  for(int i=0;i<n;i++) {
  
    int idx[11];
    double delx,dely,delz;
    const double lr = log(rho[i]);
    const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));

    get_interp_spots_d(lr,lt,ye[i],&delx,&dely,&delz,idx);

    // get prs and eps derivatives by rho and temp

    {
      int iv = 0;
      nuc_eos_C_linterp_one_d(idx,delx,dely,delz,iv,&(press[i]), &(dedrho[i]), &(dedt[i]), &(dpdrho[i]), &(dpdt[i]));
      iv = 1;
      nuc_eos_C_linterp_one_d(idx,delx,dely,delz,iv,&(eps[i]), &(dedrho[i]), &(dedt[i]), &(dpdrho[i]), &(dpdt[i]));
      iv = 7;

      nuc_eos_C_linterp_one(idx,delx,dely,delz,&(dpdeps[i]),iv);

//      dpdeps[i] = exp(dpdeps[i]);
      press[i] = exp(press[i]);
      eps2[i] = exp(eps[i]);
      eps[i] = exp(eps[i])-energy_shift;
    }

  }

  return;
}
void nuc_eos_m_kt1_dpdrhot_dedrhot(const int *restrict n_in,
				     const double *restrict rho, 
				     double *restrict temp,
				     const double *restrict ye,
				     double *restrict eps,
				     double *restrict eps2,
				     double *restrict press,
				     double *restrict dpdrho,
				     double *restrict dpdt,
				     double *restrict dedrho,
				     double *restrict dedt,
				     const double *restrict prec,
				     int *restrict keyerr,
				     int *restrict anyerr)
{

  const int n = *n_in;
  int keyerr2;

  for(int i=0;i<n;i++) {
    // check if we are fine
    keyerr2 = checkbounds(rho[i], temp[i], ye[i]);
    if(keyerr2 != 0) {
      *anyerr = 1;
      printf("Inside kt1: keyerr: %d ,%g, %g, %g",keyerr2, rho[i], temp[i], ye[i]);
    }
  }

  // Abort if there is any error in checkbounds.
  // This should never happen and the program should abort with
  // a fatal error anyway. No point in doing any further EOS calculations.
  if(*anyerr) return;

  for(int i=0;i<n;i++) {
  
    int idx[11];
    double delx,dely,delz;
    const double lr = log(rho[i]);
    const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));

    get_interp_spots_d(lr,lt,ye[i],&delx,&dely,&delz,idx);

    // get prs and eps derivatives by rho and temp

    {
      int iv = 0;
      nuc_eos_C_linterp_one_d(idx,delx,dely,delz,iv,&(press[i]), &(dedrho[i]), &(dedt[i]), &(dpdrho[i]), &(dpdt[i]));
      iv = 1;
      nuc_eos_C_linterp_one_d(idx,delx,dely,delz,iv,&(eps[i]), &(dedrho[i]), &(dedt[i]), &(dpdrho[i]), &(dpdt[i]));

      press[i] = exp(press[i]);
      eps2[i] = exp(eps[i]);
      eps[i] = exp(eps[i])-energy_shift;
    }

  }

  return;
}
void nuc_eos_m_kt1_press_eps(const int *restrict n_in,
			     const double *restrict rho, 
			     const double *restrict temp,
			     const double *restrict ye,
			     double *restrict eps,
			     double *restrict prs,
			     int *restrict keyerr,
			     int *restrict anyerr)
{

  const int n = *n_in;

  *anyerr = 0;
  for(int i=0;i<n;i++) {
    // check if we are fine
    keyerr[i] = checkbounds(rho[i], temp[i], ye[i]);
    if(keyerr[i] != 0) {
      *anyerr = 1;
      printf("Here2 keyerr = %d\n", keyerr[i]);
    }
  }

  // Abort if there is any error in checkbounds.
  // This should never happen and the program should abort with
  // a fatal error anyway. No point in doing any further EOS calculations.
  if(*anyerr) return;
  for(int i=0;i<n;i++) {
    int idx[8];
    double delx,dely,delz;
    const double xrho = log(rho[i]);
    const double xtemp = log(temp[i]);
    
    get_interp_spots(xrho,xtemp,ye[i],&delx,&dely,&delz,idx);
  
    {
      const int iv = 0;
      nuc_eos_C_linterp_one(idx,delx,dely,delz,&(prs[i]),iv);
    }

    {
      const int iv = 1;
      nuc_eos_C_linterp_one(idx,delx,dely,delz,&(eps[i]),iv);
    }
  } 

  // now get rid of ln:
  for(int i=0;i<n;i++) {
    prs[i] = exp(prs[i]);
    eps[i] = exp(eps[i]) - energy_shift;
  }


  return;
}
void nuc_eos_m_kt0_press(const int *restrict n_in,
			 const double *restrict rho, 
			 double *restrict temp,
			 const double *restrict ye,
			 const double *restrict eps,
			 double *restrict prs,
			 const double *restrict prec,
			 int *restrict keyerr,
			 int *restrict anyerr,
       int *restrict nEOScalls)
{

  const int n = *n_in;

  *anyerr = 0;
  *nEOScalls = 0;

  for(int i=0;i<n;i++) {
    // check if we are fine
    // Note that this code now requires that the
    // temperature guess be within the table bounds
    keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);
    if(keyerr[i] != 0) {
      *anyerr = 1;
    }
  }

  // Abort if there is any error in checkbounds.
  // This should never happen and the program should abort with
  // a fatal error anyway. No point in doing any further EOS calculations.
  if(*anyerr) return;

  // first must find the temperature
  for(int i=0;i<n;i++) {
    const double lr = log(rho[i]);
    const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));
    double ltout;
    const double epstot = eps[i]+energy_shift;
    if(epstot>0.0e0) {
      // this is the standard scenario; eps is larger than zero
      // and we can operate with logarithmic tables
      const double lxeps = log(epstot);
#if DEBUG
      fprintf(stderr,"%d %15.6E %15.6E %15.6E %15.6E\n",i,lr,lt,ye[i],lxeps);
      fprintf(stderr,"%d %15.6E %15.6E %15.6E %15.6E\n",i,
	      exp(lr),exp(lt),ye[i],exp(lxeps));
#endif
      int xnEOScalls = 0;
      nuc_eos_findtemp(lr,lt,ye[i],lxeps,*prec,
		       (double *restrict)(&ltout),&keyerr[i],&xnEOScalls);
      *nEOScalls += xnEOScalls;
    } else {
      // will be overwritten further down, only marks error
      keyerr[i] = 667;
    } // epstot > 0.0

    if(keyerr[i] != 0) {
      // now try negative temperature treatment
      double eps0, eps1;
      int idx[8];
      double delx,dely,delz;

      get_interp_spots_linT_low_eps(lr,temp1,ye[i],&delx,&dely,&delz,idx);
      nuc_eos_C_linterp_one_linT_low_eps(idx,delx,dely,delz,&(eps1));

      get_interp_spots_linT_low_eps(lr,temp0,ye[i],&delx,&dely,&delz,idx);
      nuc_eos_C_linterp_one_linT_low_eps(idx,delx,dely,delz,&(eps0));
      *nEOScalls += 2; 

      temp[i] = (epstot-eps0) * (temp1-temp0)/(eps1-eps0) + temp0;

      // set error codes
      *anyerr = 1;
      keyerr[i] = 668;

      // get pressure
      {
	const int iv = 0;
	get_interp_spots_linT_low(lr,temp[i],ye[i],&delx,&dely,&delz,idx);
	nuc_eos_C_linterp_one_linT_low(idx,delx,dely,delz,&(prs[i]),iv);
  *nEOScalls += 1;
      }

    } else {
      temp[i] = exp(ltout);
      int idx[8];
      double delx,dely,delz;
      get_interp_spots(lr,ltout,ye[i],&delx,&dely,&delz,idx);
      {
	const int iv = 0;
	nuc_eos_C_linterp_one(idx,delx,dely,delz,&(prs[i]),iv);
  *nEOScalls += 1;
#if DEBUG
    printf("update nuc_eos_press: nEOScalls = %i\n", *nEOScalls);
#endif  
      }
    }
  } // loop i<n

  // now get rid of ln:
  for(int i=0;i<n;i++) {
    prs[i] = exp(prs[i]);
  }

  return;
}

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 *nEOScalls) 

{

  // keyerr codes:
  // 667 -- no temperature found
  // 101 -- Y_e too high
  // 102 -- Y_e too low
  // 103 -- temp too high (if keytemp = 1)
  // 104 -- temp too low (if keytemp = 1)
  // 105 -- rho too high
  // 106 -- rho too low

  // keytemp codes:
  // 1 -- coming in with temperature
  // 0 -- coming in with eps, need to find temperature
  // 2 -- coming in with spec. entropy, need to find temperature
  //      (not currently implemented)

  int npoints = 1;
  int anyerr = 0;
  *keyerr = 0;

  // set up local vars
 
  double eps2 = 0.0;   
  double press2 = 0.0;   
  double eps = *enr;
  double lr = log(rho);
  double lt = log(*temp);
  double lenthalpy = eps;

  nuc_eos_m_kt1_press_eps(&npoints,&rho,temp,&ye,&eps2,&press2,keyerr,&anyerr);
  *nEOScalls = 1;

  double enthalpy = eps2+press2/rho;
  enthalpy = enthalpy + enthalpy*0.01;

  if(keytemp == 0) {
    double nlt = 0.0;
    int xnEOScalls = 0;
    nuc_eos_findtemp_enthalpy(lr,lt,ye,lenthalpy,rfeps,&nlt,keyerr,&xnEOScalls);
    *nEOScalls += xnEOScalls; 
#if DEBUG
    printf("update nuc_eos_press: xnEOScalls = %i\n", xnEOScalls);
#endif       
    if(*keyerr != 0) {
      printf("3 Keyerr: %d, Temp: %g Rho: %g Eps: %g Ye: %g\n",*keyerr,*temp,rho,eps,ye); 
    }
    lt = nlt;

    *temp = exp(lt);
  } else if(keytemp == 1) {
    
  }

  int idx[8];
  double delx,dely,delz,lpress,leps;

  get_interp_spots(lr,lt,ye,&delx,&dely,&delz,idx);
  nuc_eos_C_linterp_one(idx,delx,dely,delz,&lpress,0);  
  nuc_eos_C_linterp_one(idx,delx,dely,delz,&leps,1);  
  *nEOScalls += 1;
  // assign results

  *enr2 = exp(leps) - energy_shift;
  *press= exp(lpress);
  return;
}