/*@@
   @file      EOS_GeneralHybrid.F90
   @date      
   @author    Ian Hawke, Christian D. Ott
   @desc 
      Routines to calculate the EOS used by Dimmelmeier et al. 
      in supernova core collapse simulations.
   @enddesc 
 @@*/

#include "cctk.h"
#include "cctk_Parameters.h"

subroutine EOS_GeneralHybrid_Pressure(nelems, rho, eps, press)

  USE EOS_GP_Scalars   
  USE EOS_GeneralHybrid_Scalars

  implicit none

  DECLARE_CCTK_PARAMETERS

  CCTK_INT, intent(in) :: nelems
  CCTK_REAL, dimension(nelems), intent(in) :: rho
  CCTK_REAL, dimension(nelems), intent(in) :: eps
  CCTK_REAL, dimension(nelems), intent(out) :: press

  CCTK_REAL, allocatable, dimension(:) :: p_poly
  CCTK_REAL, allocatable, dimension(:) :: p_th
  CCTK_REAL, allocatable, dimension(:) :: local_eos_gamma
  CCTK_REAL, allocatable, dimension(:) :: local_eos_k_cgs

  CCTK_REAL zero

  allocate(p_poly(1:nelems))
  allocate(p_th(1:nelems))
  allocate(local_eos_gamma(1:nelems))
  allocate(local_eos_k_cgs(1:nelems))

  zero = 0.d0

  where (rho > rho_nuc) 
     local_eos_gamma = eos_gamma_supernuclear
     local_eos_k_cgs = eos_k_supernuclear_cgs
  elsewhere
     local_eos_gamma = eos_gamma
     local_eos_k_cgs = eos_k_cgs
  end where

  p_poly = p_geom_factor * local_eos_k_cgs *  &
   (rho * rho_geom_factor_inv)**local_eos_gamma

  p_th = - p_geom_factor * local_eos_k_cgs * (eos_gamma_th - 1.d0) /  &
       (local_eos_gamma - 1.d0) * (rho * rho_geom_factor_inv)**local_eos_gamma + &
       (eos_gamma_th - 1.d0) * rho * eps - &
       (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / &
       (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * &
       p_geom_factor * eos_k_cgs * rho_geom_factor_inv**eos_gamma * &
       rho_nuc**(eos_gamma - 1.d0) * rho
  
  p_th = max(zero, p_th)

  press = p_poly + p_th

  deallocate(p_poly)
  deallocate(p_th)
  deallocate(local_eos_gamma)
  deallocate(local_eos_k_cgs)


end subroutine EOS_GeneralHybrid_Pressure

subroutine EOS_GeneralHybrid_DPDRho(nelems, rho, eps, dpdrho)

  USE EOS_GP_Scalars
  USE EOS_GeneralHybrid_Scalars

  implicit none
  DECLARE_CCTK_PARAMETERS

  CCTK_INT, intent(in) :: nelems
  CCTK_REAL, dimension(nelems), intent(in) :: rho
  CCTK_REAL, dimension(nelems), intent(in) :: eps
  CCTK_REAL, dimension(nelems), intent(out) :: dpdrho

  CCTK_REAL, allocatable,dimension(:) :: d_p_poly
  CCTK_REAL, allocatable,dimension(:) :: d_p_th_1
  CCTK_REAL, allocatable,dimension(:) :: d_p_th_2
  CCTK_REAL, allocatable,dimension(:) :: local_eos_gamma
  CCTK_REAL, allocatable,dimension(:) :: local_eos_k_cgs

  CCTK_REAL zero

  allocate(d_p_poly(1:nelems))
  allocate(d_p_th_1(1:nelems))
  allocate(d_p_th_2(1:nelems))
  allocate(local_eos_gamma(1:nelems))
  allocate(local_eos_k_cgs(1:nelems))

  
  zero = 0.d0

  where (rho > rho_nuc) 
     local_eos_gamma = eos_gamma_supernuclear
     local_eos_k_cgs = eos_k_supernuclear_cgs
  elsewhere
     local_eos_gamma = eos_gamma
     local_eos_k_cgs = eos_k_cgs
  end where

  d_p_poly = local_eos_gamma * p_geom_factor * local_eos_k_cgs * &
       rho**(local_eos_gamma - 1.d0) * rho_geom_factor_inv**local_eos_gamma

  d_p_th_1 = - local_eos_gamma * p_geom_factor * local_eos_k_cgs * &
       (eos_gamma_th - 1.d0) / (local_eos_gamma - 1.d0) * &
       rho**(local_eos_gamma - 1.d0) * rho_geom_factor_inv**local_eos_gamma

  d_p_th_2 = (eos_gamma_th - 1.d0) * eps &
       - (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / &
       (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * &
       p_geom_factor * eos_k_cgs * rho_geom_factor_inv**eos_gamma * &
       rho_nuc**(eos_gamma - 1.d0)

!  d_p_th_1 = max(d_p_th_1, zero)
!  d_p_th_2 = max(d_p_th_2, zero)
      
  dpdrho = d_p_poly + d_p_th_1 + d_p_th_2

  deallocate(d_p_poly)
  deallocate(d_p_th_1)
  deallocate(d_p_th_2)
  deallocate(local_eos_gamma)
  deallocate(local_eos_k_cgs)


end subroutine EOS_GeneralHybrid_DPDRho


subroutine EOS_GeneralHybrid_cs2(nelems, rho, eps, cs2)

! This is an ugly and quick hack.
! Uses way to many operations and memory.

  USE EOS_GP_Scalars
  USE EOS_GeneralHybrid_Scalars

  implicit none
  DECLARE_CCTK_PARAMETERS

  CCTK_INT, intent(in) :: nelems
  CCTK_REAL, dimension(nelems), intent(in) :: rho
  CCTK_REAL, dimension(nelems), intent(in) :: eps
  CCTK_REAL, dimension(nelems), intent(out) :: cs2
  CCTK_REAL, allocatable,dimension(:) :: p_poly
  CCTK_REAL, allocatable,dimension(:) :: p_th
  CCTK_REAL, allocatable,dimension(:) :: pressure

  CCTK_REAL, allocatable,dimension(:) :: local_eos_gamma
  CCTK_REAL, allocatable,dimension(:) :: local_eos_k_cgs

  CCTK_REAL zero

  allocate(p_poly(1:nelems))
  allocate(p_th(1:nelems))
  allocate(pressure(1:nelems))
  allocate(local_eos_gamma(1:nelems))
  allocate(local_eos_k_cgs(1:nelems))


  zero = 0.d0

  where (rho > rho_nuc) 
     local_eos_gamma = eos_gamma_supernuclear
     local_eos_k_cgs = eos_k_supernuclear_cgs
  elsewhere
     local_eos_gamma = eos_gamma
     local_eos_k_cgs = eos_k_cgs
  end where

! First calculate the pressure 
  p_poly = p_geom_factor * local_eos_k_cgs *  &
   (rho * rho_geom_factor_inv)**local_eos_gamma

  p_th = - p_geom_factor * local_eos_k_cgs * (eos_gamma_th - 1.d0) /  &
       (local_eos_gamma - 1.d0) * (rho * rho_geom_factor_inv)**local_eos_gamma + &
       (eos_gamma_th - 1.d0) * rho * eps - &
       (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / &
       (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * &
       p_geom_factor * eos_k_cgs * rho_geom_factor_inv**eos_gamma * &
       rho_nuc**(eos_gamma - 1.d0) * rho

  pressure = p_poly + p_th
  
  p_th = max(zero, p_th)

! This may look incorrect; It's in Harry's code; in his thesis and I worked it out as well.
! Should be okay...
  cs2 = (local_eos_gamma * p_poly + eos_gamma_th * p_th) / &
       rho / (1.0d0 + eps + pressure/rho)

  deallocate(pressure)
  deallocate(p_poly)
  deallocate(p_th)
  deallocate(local_eos_gamma)
  deallocate(local_eos_k_cgs)

end subroutine EOS_GeneralHybrid_cs2



subroutine EOS_GeneralHybrid_DPDIE(nelems, rho, eps, dpdie)


  USE EOS_GP_Scalars
  USE EOS_GeneralHybrid_Scalars

  implicit none
  DECLARE_CCTK_PARAMETERS

  CCTK_INT, intent(in) :: nelems
  CCTK_REAL, dimension(nelems), intent(in) :: rho
  CCTK_REAL, dimension(nelems), intent(in) :: eps
  CCTK_REAL, dimension(nelems), intent(out) :: dpdie

  dpdie = (eos_gamma_th - 1.d0) * rho

end subroutine EOS_GeneralHybrid_DPDIE