* update nuc_eos with some changes in the low-density regime.

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

1 files changed, 32 insertions, 13 deletions

diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90
index b576ce0..64a78a5 100644
--- a/src/nuc_eos/nuc_eos.F90
+++ b/src/nuc_eos/nuc_eos.F90
@@ -160,6 +160,25 @@ subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
 
 end subroutine nuc_eos_full
 
+subroutine nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+
+  implicit none
+  real*8 xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe
+  real*8,parameter:: polyK1 = 5.0d14
+  real*8,parameter :: polyGamma = 1.33333333333333d0
+  integer keytemp
+
+  xprs=polyK1*xrho**(polyGamma)
+  xenr=xprs/xrho/(polyGamma-1.d0)
+  xcs2=polyGamma*xprs/xrho
+  xdpderho=0.0d0
+  xdpdrhoe=0.0d0
+
+  xenr=max(8.0d18,xenr)
+  xenr=min(3.0d20,xenr)
+
+end subroutine nuc_poly_eos
+
 subroutine nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
 
   implicit none
@@ -212,9 +231,9 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
      call CCTK_WARN (0, "nuc_eos: rho > rhomax")
   endif
 
-  if(xrho.lt.eos_rhomin*1.2d0) then
-     call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
-     xent = 4.0d0
+  if(xrho.lt.eos_rhomin*1.126d0) then
+     call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+     xent = 0.0d0
      return
   endif
 
@@ -234,8 +253,8 @@ subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
      endif
      
      if(xtemp.lt.eos_tempmin) then
-        call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
-        xent = 4.0d0
+        call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+        xent = 0.0d0
         return
      endif
   endif
@@ -335,8 +354,8 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
      call CCTK_WARN (0, "nuc_eos: rho > rhomax")
   endif
 
-  if(xrho.lt.eos_rhomin*1.2d0) then
-     call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+  if(xrho.lt.eos_rhomin*1.126d0) then
+     call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
      return
   endif
 
@@ -354,9 +373,9 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
      if(xtemp.gt.eos_tempmax) then
         call CCTK_WARN (0, "nuc_eos: temp > tempmax")
      endif
-     
+
      if(xtemp.lt.eos_tempmin) then
-        call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+        call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
         return
      endif
   endif
@@ -414,8 +433,8 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
      call CCTK_WARN (0, "nuc_eos: rho > rhomax")
   endif
 
-  if(xrho.lt.eos_rhomin*1.2d0) then
-     call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+  if(xrho.lt.eos_rhomin*1.126d0) then
+     call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
      return
   endif
 
@@ -433,9 +452,9 @@ subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
      if(xtemp.gt.eos_tempmax) then
         call CCTK_WARN (0, "nuc_eos: temp > tempmax")
      endif
-     
+
      if(xtemp.lt.eos_tempmin) then
-        call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+        call nuc_poly_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
         return
      endif
   endif