* add aliased routine to compute deps/drho|p and deps/dp|rho

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@57 8e189c6b-2ab8-4400-aa02-70a9cfce18b9

2 files changed, 116 insertions, 0 deletions

diff --git a/interface.ccl b/interface.ccl
index 55d6d58..57ef222 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -138,8 +138,23 @@ void FUNCTION EOS_Omni_short(CCTK_INT IN eoskey,                \
 
 PROVIDES FUNCTION EOS_Omni_short WITH EOS_Omni_EOS_short LANGUAGE Fortran
 
+################################################################################
+# the following routines are needed for MHD con2prim using the 2D Z-P scheme
 
+void FUNCTION EOS_Omni_DEpsByDRho_DEpsByDPress(CCTK_INT IN eoskey,         \
+			     CCTK_INT IN havetemp,              \
+			     CCTK_REAL IN rf_precision,         \
+                             CCTK_INT IN npoints,               \
+			     CCTK_REAL IN ARRAY rho,     	\
+			     CCTK_REAL INOUT ARRAY eps,  	\
+			     CCTK_REAL INOUT ARRAY temp, 	\
+			     CCTK_REAL IN ARRAY ye,      	\
+			     CCTK_REAL OUT ARRAY DEpsByDRho,    \
+			     CCTK_REAL OUT ARRAY DEpsByDPress,  \
+			     CCTK_INT OUT ARRAY keyerr,   	\
+                             CCTK_INT OUT anyerr)
 
+PROVIDES FUNCTION EOS_Omni_DEpsByDRho_DEpsByDPress WITH EOS_Omni_EOS_DEpsByDRho_DEpsByDPress LANGUAGE Fortran
 ################################################################################
 # the following routines are needed for MHD con2prim based on the Harm Code
 
diff --git a/src/EOS_Omni_MultiVarCalls.F90 b/src/EOS_Omni_MultiVarCalls.F90
index 664a31f..a989fd3 100644
--- a/src/EOS_Omni_MultiVarCalls.F90
+++ b/src/EOS_Omni_MultiVarCalls.F90
@@ -202,3 +202,104 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
    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_cgs
+  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) = press_gf * poly_k_cgs * &
+                   (rho(i)*inv_rho_gf)**(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) * press_gf * poly_k_cgs *  &
+                poly_gamma * inv_rho_gf *        &
+                (rho(i)*inv_rho_gf) ** (poly_gamma - 1.d0)
+        enddo
+  case (2)
+     ! gamma-law EOS                                                                                       
+        if(keytemp.eq.1) then
+           do i=1,npoints
+              eps(i) = press_gf * gl_k_cgs * &
+                   (rho(i)*inv_rho_gf)**(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_cgs = hybrid_k2_cgs
+           else
+              hybrid_local_gamma = hybrid_gamma1
+              hybrid_local_k_cgs = hybrid_k1_cgs
+           endif
+           hybrid_dp_poly = hybrid_local_gamma * press_gf *                 &
+                hybrid_local_k_cgs * rho(i)**(hybrid_local_gamma - 1.0d0) * &
+                inv_rho_gf**hybrid_local_gamma
+
+           hybrid_dp_th1 = - hybrid_local_gamma * press_gf * hybrid_local_k_cgs * &
+                (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) *      &
+                rho(i)**(hybrid_local_gamma - 1.d0) * inv_rho_gf**hybrid_local_gamma
+
+           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) *                    &
+                press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 *               &
+                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