#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"

! eoskey:
! 1 --- polytropic EOS
! 2 --- gamma-law EOS
! 3 --- hybrid EOS
! 4 --- finite-T microphysical NSE EOS


subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
     rho,eps,temp,ye,press,entropy,cs2,dedt,dpderho,dpdrhoe,munu,&
     keyerr,anyerr)

  use EOS_Omni_Module
  implicit none
  DECLARE_CCTK_PARAMETERS

  CCTK_INT, intent(in)     :: eoskey,keytemp,npoints
  CCTK_INT, intent(out)    :: keyerr(npoints)
  CCTK_INT, intent(out)    :: anyerr
  CCTK_REAL, intent(in)    :: rf_precision
  CCTK_REAL, intent(in)    :: rho(npoints),ye(npoints)
  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
  CCTK_REAL, intent(out)   :: press(npoints)
  CCTK_REAL, intent(inout) :: entropy(npoints)
  CCTK_REAL, intent(out)   :: cs2(npoints)
  CCTK_REAL, intent(out)   :: dedt(npoints)
  CCTK_REAL, intent(out)   :: dpderho(npoints)
  CCTK_REAL, intent(out)   :: dpdrhoe(npoints)
  CCTK_REAL, intent(out)   :: munu(npoints)

  ! local vars
  integer          :: i
  character(256)   :: warnstring
  ! temporary vars for nuc_eos
  real*8           :: xrho,xye,xtemp,xenr,xent
  real*8           :: xprs,xmunu,xcs2
  real*8           :: xdedt,xdpderho,xdpdrhoe

  if(eoskey.ne.4) then
     call CCTK_WARN(0,"EOS_Omni_EOS_short currently does not work for this eoskey")
  endif

  anyerr    = 0
  keyerr(:) = 0

  do i=1,npoints

     xrho = rho(i) * inv_rho_gf
     xtemp = temp(i)
     xye = ye(i)
     xenr = eps(i) * inv_eps_gf
     xent = entropy(i)
     call nuc_eos_short(xrho,xtemp,xye,xenr,xprs,& 
          xent,xcs2,xdedt,xdpderho,xdpdrhoe,xmunu,&
          keytemp,keyerr(i),rf_precision)
     
     if(keyerr(i).ne.0) then
        anyerr = 1
     endif
     
     if(keytemp.eq.1) then
        eps(i) = xenr * eps_gf
     else if(keytemp.eq.2) then
        eps(i) = xenr * eps_gf
        temp(i) = xtemp
     else
        temp(i) = xtemp
     endif
     
     press(i)     = xprs * press_gf
     entropy(i)   = xent
     cs2(i)       = xcs2
     dedt(i)      = xdedt * eps_gf
     dpderho(i)   = xdpderho * press_gf * inv_eps_gf
     dpdrhoe(i)   = xdpdrhoe * press_gf * inv_rho_gf
     munu(i)      = xmunu
     
  enddo

end subroutine EOS_Omni_EOS_short


subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
     rho,eps,temp,ye,dpderho,dpdrhoe,keyerr,anyerr)

  use EOS_Omni_Module
  implicit none
  DECLARE_CCTK_PARAMETERS

  CCTK_INT, intent(in)     :: eoskey,keytemp,npoints
  CCTK_INT, intent(out)    :: keyerr(npoints)
  CCTK_INT, intent(out)    :: anyerr
  CCTK_REAL, intent(in)    :: rf_precision
  CCTK_REAL, intent(in)    :: rho(npoints),ye(npoints)
  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
  CCTK_REAL, intent(out)   :: dpderho(npoints)
  CCTK_REAL, intent(out)   :: dpdrhoe(npoints)

  ! local vars
  integer          :: i
  character(256)   :: warnstring
  real*8           :: hybrid_local_gamma
  real*8           :: hybrid_local_k
  real*8           :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2
  ! temporary vars for nuc_eos
  real*8           :: xrho,xye,xtemp,xenr,xent
  real*8           :: xprs,xmunu,xcs2
  real*8           :: xdedt,xdpderho,xdpdrhoe

  anyerr    = 0
  keyerr(:) = 0

  select case (eoskey)
  case (1)
     ! polytropic EOS                                                                                      
        if(keytemp.eq.1) then
           do i=1,npoints
              eps(i) = poly_k * &
                   rho(i)**(poly_gamma) / &
                   (poly_gamma - 1.0d0) / rho(i)
           enddo
        endif
        do i=1,npoints
           dpdrhoe(i) = poly_k *  &
                poly_gamma *        &
                rho(i) ** (poly_gamma - 1.d0)
           dpderho(i) = 0.0d0
        enddo
  case (2)
     ! gamma-law EOS                                                                                       
        if(keytemp.eq.1) then
           do i=1,npoints
              eps(i) = gl_k * &
                   rho(i)**(gl_gamma) / &
                   (gl_gamma - 1.0d0) / rho(i)
           enddo
        endif
        do i=1,npoints
           dpdrhoe(i) = (gl_gamma-1.0d0) * &
                eps(i)
           dpderho(i) = (gl_gamma - 1.0d0) * &
                rho(i)
        enddo
  case (3)
     ! hybrid EOS                                                                                          
        do i=1,npoints
           if(rho(i).gt.hybrid_rho_nuc) then
              hybrid_local_gamma = hybrid_gamma2
              hybrid_local_k = hybrid_k2
           else
              hybrid_local_gamma = hybrid_gamma1
              hybrid_local_k = hybrid_k1
           endif
           hybrid_dp_poly = hybrid_local_gamma *                 &
                hybrid_local_k * rho(i)**(hybrid_local_gamma - 1.0d0)
                

           hybrid_dp_th1 = - hybrid_local_gamma * hybrid_local_k * &
                (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) *      &
                rho(i)**(hybrid_local_gamma - 1.d0)

           hybrid_dp_th2  = (hybrid_gamma_th - 1.d0) * eps(i)                        &
                - (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) /  &
                (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) *                    &
                hybrid_k1 *                &
                hybrid_rho_nuc**(hybrid_gamma1 - 1.d0)

           dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2
           dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
        enddo
    case (4)
        do i=1,npoints
           xrho = rho(i) * inv_rho_gf
           xtemp = temp(i)
           xye = ye(i)
           xenr = eps(i) * inv_eps_gf
           call nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr, &
                xdpderho, xdpdrhoe,&
                keytemp,keyerr(i),rf_precision)

           if(keyerr(i).ne.0) then
              anyerr = 1
           endif

           if(keytemp.eq.1) then
              eps(i) = xenr * eps_gf
           else
              temp(i) = xtemp
           endif

           dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf
           dpderho(i) = xdpderho * press_gf * inv_eps_gf
        enddo
   case DEFAULT
      write(warnstring,*) "eoskey ",eoskey," not implemented!"
      call CCTK_WARN(0,warnstring)
   end select

 end subroutine EOS_Omni_EOS_dpderho_dpdrhoe


 subroutine EOS_Omni_EOS_DEpsByDRho_DEpsByDPress(eoskey,keytemp,rf_precision,npoints,&
     rho,eps,temp,ye,depsdrho,depsdpress,keyerr,anyerr)

  use EOS_Omni_Module
  implicit none
  DECLARE_CCTK_PARAMETERS

  CCTK_INT, intent(in)     :: eoskey,keytemp,npoints
  CCTK_INT, intent(out)    :: keyerr(npoints)
  CCTK_INT, intent(out)    :: anyerr
  CCTK_REAL, intent(in)    :: rf_precision
  CCTK_REAL, intent(in)    :: rho(npoints),ye(npoints)
  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
  CCTK_REAL, intent(out)   :: depsdrho(npoints)
  CCTK_REAL, intent(out)   :: depsdpress(npoints)

  ! local vars
  integer          :: i
  character(256)   :: warnstring
  real*8           :: hybrid_local_gamma
  real*8           :: hybrid_local_k
  real*8           :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2
  ! temporary vars for nuc_eos
  real*8           :: xrho,xye,xtemp,xenr,xent
  real*8           :: xprs,xmunu,xcs2
  real*8           :: xdedt,xdpderho,xdpdrhoe


  anyerr    = 0
  keyerr(:) = 0

  select case (eoskey)
  case (1)
     ! polytropic EOS                                                                                      
        if(keytemp.eq.1) then
           do i=1,npoints
              eps(i) = poly_k * &
                   rho(i)**(poly_gamma) / &
                   (poly_gamma - 1.0d0) / rho(i)
           enddo
        endif
        do i=1,npoints
           depsdpress(i) = 1.0d0/(poly_gamma - 1.0d0)/rho(i)
           depsdrho(i) = depsdpress(i) * poly_k *  &
                poly_gamma *         &
                rho(i) ** (poly_gamma - 1.d0)
        enddo
  case (2)
     ! gamma-law EOS                                                                                       
        if(keytemp.eq.1) then
           do i=1,npoints
              eps(i) = gl_k * &
                   rho(i)**(gl_gamma) / &
                   (gl_gamma - 1.0d0) / rho(i)
           enddo
        endif
        do i=1,npoints
           depsdpress(i) = 1.0/( (gl_gamma - 1.0d0) * &
                rho(i))
           depsdrho(i) = -eps(i)/rho(i)
       enddo
  case (3)
     ! hybrid EOS                                                                                          
        do i=1,npoints
           if(rho(i).gt.hybrid_rho_nuc) then
              hybrid_local_gamma = hybrid_gamma2
              hybrid_local_k = hybrid_k2
           else
              hybrid_local_gamma = hybrid_gamma1
              hybrid_local_k = hybrid_k1
           endif
           hybrid_dp_poly = hybrid_local_gamma *                 &
                hybrid_local_k * rho(i)**(hybrid_local_gamma - 1.0d0)
                

           hybrid_dp_th1 = - hybrid_local_gamma * hybrid_local_k * &
                (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) *      &
                rho(i)**(hybrid_local_gamma - 1.d0)

           hybrid_dp_th2  = (hybrid_gamma_th - 1.d0) * eps(i)                        &
                - (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) /  &
                (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) *                    &
                hybrid_k1 *               &
                hybrid_rho_nuc**(hybrid_gamma1 - 1.d0)

           xdpdrhoe = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2
           xdpderho = (hybrid_gamma_th - 1.0d0) * rho(i)
           depsdpress(i) = 1.0 / xdpderho
           depsdrho(i) = - xdpdrhoe * depsdpress(i)
        enddo
    case (4)
      write(warnstring,*) "depsdrho and depsdpress not implemented yet for hot nuclear EOS"
      call CCTK_WARN(0,warnstring)
   case DEFAULT
      write(warnstring,*) "eoskey ",eoskey," not implemented!"
      call CCTK_WARN(0,warnstring)
   end select

 end subroutine EOS_Omni_EOS_DEpsByDRho_DEpsByDPress