/*@@
   @file      EOS_GP.F90
   @date      Mon Mar 14 16:42:15 2005
   @author    Ian Hawke
   @desc 
   The functions that actually set the polytropic EOS
   @enddesc 
 @@*/

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

subroutine EOS_GP_Pressure(nelems, rho, press)

  USE EOS_GP_Scalars

  implicit none

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

  press = p_geom_factor * eos_k_cgs * &
       (rho * rho_geom_factor_inv) ** eos_gamma_local  

end subroutine EOS_GP_Pressure

subroutine EOS_GP_IntEn(nelems, rho, IntEn)

  USE EOS_GP_Scalars

  implicit none

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

  IntEn = p_geom_factor * eos_k_cgs * &
       (rho * rho_geom_factor_inv) ** eos_gamma_local / &
       (rho * (eos_gamma_local - 1.d0) )

end subroutine EOS_GP_IntEn

subroutine EOS_GP_DPDRho(nelems, rho, dpdrho)

  USE EOS_GP_Scalars

  implicit none

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

  dpdrho = p_geom_factor * eos_k_cgs * &
           eos_gamma_local * rho_geom_factor_inv * &
       (rho * rho_geom_factor_inv) ** (eos_gamma_local - 1.d0)

end subroutine EOS_GP_DPDRho

subroutine EOS_GP_DPDIE(nelems, rho, dpdie)

  USE EOS_GP_Scalars

  implicit none

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

  dpdie = 0.d0

end subroutine EOS_GP_DPDIE

subroutine EOS_GP_cs2(nelems, rho, cs2)

  USE EOS_GP_Scalars

  implicit none

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

  cs2 = (p_geom_factor * eos_k_cgs * eos_gamma_local * rho_geom_factor_inv * &
         (rho * rho_geom_factor_inv) ** (eos_gamma_local - 1.d0)) / &
       (1.d0 + p_geom_factor * eos_k_cgs * eos_gamma_local * &
       rho_geom_factor_inv / (eos_gamma_local - 1.d0) * &
       (rho * rho_geom_factor_inv) ** (eos_gamma_local - 1.d0))
  
end subroutine EOS_GP_cs2

 /*@@
   @routine    Inverse functions
   @date       Mon Apr  4 11:03:04 2005
   @author     Ian Hawke
   @desc 
   Give P find the results...
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


subroutine EOS_GP_Inv_Rho(nelems, press, rho)

  USE EOS_GP_Scalars

  implicit none

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

  rho = rho_geom_factor * &
       (press / p_geom_factor / eos_k_cgs) ** (1.d0 / eos_gamma_local)

end subroutine EOS_GP_Inv_Rho

subroutine EOS_GP_Inv_IntEn(nelems, press, IntEn)

  USE EOS_GP_Scalars

  implicit none

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

  IntEn = press / &
       ( (eos_gamma_local - 1.d0) * rho_geom_factor * &
       (press / p_geom_factor / eos_k_cgs) ** (1.d0 / eos_gamma_local) )

end subroutine EOS_GP_Inv_IntEn

subroutine EOS_GP_Inv_DPDRho(nelems, press, dpdrho)

  USE EOS_GP_Scalars

  implicit none

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

  dpdrho = p_geom_factor * eos_k_cgs * &
           eos_gamma_local * rho_geom_factor_inv * &
       (press / (p_geom_factor * eos_k_cgs) ) ** &
       ( (eos_gamma_local - 1.d0) / eos_gamma_local )

end subroutine EOS_GP_Inv_DPDRho

subroutine EOS_GP_Inv_DPDIE(nelems, press, dpdie)

  USE EOS_GP_Scalars

  implicit none

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

  dpdie = 0.d0

end subroutine EOS_GP_Inv_DPDIE

subroutine EOS_GP_Inv_cs2(nelems, press, cs2)

  USE EOS_GP_Scalars

  implicit none

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

  cs2 = (p_geom_factor * eos_k_cgs * eos_gamma_local * rho_geom_factor_inv * &
         (press / (p_geom_factor * eos_k_cgs) ) ** &
         ( (eos_gamma_local - 1.d0) / eos_gamma_local ) ) / &
       (1.d0 + p_geom_factor * eos_k_cgs * eos_gamma_local * &
       rho_geom_factor_inv / (eos_gamma_local - 1.d0) * &
       (press / (p_geom_factor * eos_k_cgs) ) ** &
       ( (eos_gamma_local - 1.d0) / eos_gamma_local ) )
  
end subroutine EOS_GP_Inv_cs2