thorns for general eos interface (moved from Whisky_Dev repositories)

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOSG_Polytrope/trunk@2 a192deed-9f5d-41fb-bb60-7f59df5b1b8d

1 files changed, 189 insertions, 0 deletions

diff --git a/src/EOS_GP.F90 b/src/EOS_GP.F90
new file mode 100755
index 0000000..283b9ac
--- /dev/null
+++ b/src/EOS_GP.F90
@@ -0,0 +1,189 @@
+ /*@@
+   @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
+