[Commits] [svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos/ (Rev. 37)

[cott at tapir.caltech.edu]

User: cott
Date: 2011/02/06 11:36 PM

Modified:
 /trunk/src/nuc_eos/
  nuc_eos.F90

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

File Changes:

Directory: /trunk/src/nuc_eos/
==============================

File [modified]: nuc_eos.F90
Delta lines: +32 -13
===================================================================
--- trunk/src/nuc_eos/nuc_eos.F90	2011-02-05 17:07:51 UTC (rev 36)
+++ trunk/src/nuc_eos/nuc_eos.F90	2011-02-07 05:36:36 UTC (rev 37)
@@ -160,6 +160,25 @@
 
 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 @@
      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 @@
      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 @@
      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 @@
      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 @@
      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 @@
      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