* add finite-T EOS routines (not called yet; but code compiles!)

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

8 files changed, 1043 insertions, 0 deletions

diff --git a/src/nuc_eos/bisection.F90 b/src/nuc_eos/bisection.F90
new file mode 100644
index 0000000..b8edec2
--- /dev/null
+++ b/src/nuc_eos/bisection.F90
@@ -0,0 +1,82 @@
+subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect)
+
+  use eosmodule
+
+  implicit none
+
+  real*8 lr,lt0,y,eps0,lt
+  integer keyerrt
+
+  integer keybisect
+  ! 1 -> operate on energy
+  ! 2 -> operate on entropy
+
+  !temporary vars
+  real*8 lt1,lt2,ltmin,ltmax
+  real*8 f1,f2,fmid,dlt,ltmid
+  real*8 d1,d2,d3,tol
+  real*8 f1a,f2a
+
+  integer bcount,i,itmax,maxbcount
+
+  real*8 bivar
+
+  bcount = 0
+  maxbcount = 80
+
+  tol=1.d-9 ! need to find energy to less than 1 in 10^-9
+  itmax=50
+
+  ltmax=logtemp(ntemp)
+  ltmin=logtemp(1)
+
+  lt = lt0
+  lt1 = dlog10(min(10.0d0**ltmax,1.10d0*(10.0d0**lt0)))
+  lt2 = dlog10(max(10.0d0**ltmin,0.90d0*(10.0d0**lt0)))
+
+  call findthis(lr,lt1,y,f1a,bivar,d1,d2,d3)
+  call findthis(lr,lt2,y,f2a,bivar,d1,d2,d3)
+
+  f1=f1a-eps0
+  f2=f2a-eps0
+
+  keyerrt=0
+  
+
+  do while(f1*f2.ge.0.0d0)
+     bcount=bcount+1
+     lt1=dlog10(min(10.0d0**ltmax,1.2d0*(10.0d0**lt1)))
+     lt2=dlog10(max(10.0d0**ltmin,0.8d0*(10.0d0**lt2)))
+     call findthis(lr,lt1,y,f1a,bivar,d1,d2,d3)
+     call findthis(lr,lt2,y,f2a,bivar,d1,d2,d3)
+     f1=f1a-eps0
+     f2=f2a-eps0
+     if(bcount.ge.maxbcount) then
+        keyerrt=667
+        return
+     endif
+  enddo
+
+  if(f1.lt.0.0d0) then
+     lt=lt1
+     dlt=dlog10((10.0D0**lt2)-(10.0d0**lt1))
+  else
+     lt=lt2
+     dlt=dlog10((10.0D0**lt1)-(10.0d0**lt2))
+  endif
+
+  do i=1,itmax
+     dlt=dlog10((10.0d0**dlt)*0.5d0)
+     ltmid=dlog10(10.0d0**lt+10.0d0**dlt)
+     call findthis(lr,ltmid,y,f2a,bivar,d1,d2,d3)
+     fmid=f2a-eps0
+     if(fmid.le.0.0d0) lt=ltmid
+     if(abs(1.0d0-f2a/eps0).lt.tol) then
+        lt=ltmid
+        return
+     endif
+  enddo
+
+
+
+end subroutine bisection
diff --git a/src/nuc_eos/eosmodule.F90 b/src/nuc_eos/eosmodule.F90
new file mode 100644
index 0000000..ed5e7c6
--- /dev/null
+++ b/src/nuc_eos/eosmodule.F90
@@ -0,0 +1,60 @@
+ module eosmodule
+
+   implicit none
+   
+   integer,save :: nrho,ntemp,nye
+
+   integer,save :: warn_from !warn from given reflevel
+
+   real*8 :: energy_shift = 0.0d0
+   
+   real*8 :: precision = 1.0d-9
+
+! min-max values:
+   real*8 :: eos_rhomin,eos_rhomax
+   real*8 :: eos_yemin,eos_yemax
+   real*8 :: eos_tempmin,eos_tempmax
+
+   real*8 :: t_max_hack = 240.0d0
+
+! basics
+   integer, parameter :: nvars = 19
+   real*8,allocatable :: alltables(:,:,:,:)
+  ! index variable mapping:
+  !  1 -> logpress
+  !  2 -> logenergy
+  !  3 -> entropy
+  !  4 -> munu
+  !  5 -> cs2
+  !  6 -> dedT
+  !  7 -> dpdrhoe
+  !  8 -> dpderho
+  !  9 -> muhat
+  ! 10 -> mu_e
+  ! 11 -> mu_p
+  ! 12 -> mu_n
+  ! 13 -> xa
+  ! 14 -> xh
+  ! 15 -> xn
+  ! 16 -> xp
+  ! 17 -> abar
+  ! 18 -> zbar
+  ! 19 -> gamma
+
+   real*8,allocatable,save :: logrho(:)
+   real*8,allocatable,save :: logtemp(:)
+   real*8,allocatable,save :: ye(:)
+
+! constants
+   real*8,save :: mev_to_erg = 1.60217733d-6
+   real*8,save :: amu_cgs = 1.66053873d-24
+   real*8,save :: amu_mev = 931.49432d0
+   real*8,save :: pi = 3.14159265358979d0
+   real*8,save :: ggrav = 6.672d-8
+   real*8,save :: temp_mev_to_kelvin = 1.1604447522806d10
+   real*8,save :: clight = 2.99792458d10
+   real*8,save :: kb_erg = 1.380658d-16
+   real*8,save :: kb_mev = 8.61738568d-11   
+
+
+ end module eosmodule
diff --git a/src/nuc_eos/findtemp.F90 b/src/nuc_eos/findtemp.F90
new file mode 100644
index 0000000..d208a5f
--- /dev/null
+++ b/src/nuc_eos/findtemp.F90
@@ -0,0 +1,207 @@
+subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
+
+  use eosmodule
+
+  implicit none
+
+  real*8 lr,lt0,y,epsin
+  real*8 eps,lt,ldt
+  real*8 tol
+  real*8 d1,d2,d3
+  real*8 eps0,eps1,lt1
+
+  real*8 ltn,ltmax,ltmin
+  real*8 tinput,rfeps
+
+  integer :: rl = 0
+  integer itmax,i,keyerrt
+  integer ii,jj,kk
+  
+  keyerrt=0
+
+  tol=rfeps ! need to find energy to less than 1 in 10^-10
+  itmax=20 ! use at most 20 iterations, then bomb
+
+  lt=lt0
+  lt1=lt 
+
+  eps0=epsin
+  eps1=eps0
+  
+
+  ltmax=logtemp(ntemp)
+  ltmin=logtemp(1)
+
+  ! Note: We are using Ewald's Lagrangian interpolator here!
+
+  !preconditioning 1: do we already have the right temperature?
+  call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3)
+  if (abs(eps-eps0).lt.tol*abs(eps0)) then
+     return
+  endif
+  lt1=lt
+  eps1=eps
+ 
+!  write(*,"(i4,1P12E19.10)") 0,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0,d2
+!  write(*,"(i4,1P12E19.10)") 0,lr,lt0,y,ltmin,ltmax
+
+  do i=1,itmax
+     !d2 is the derivative deps/dlogtemp;
+     ldt = -(eps - eps0)/d2 
+!     if(ldt.gt.0.d0) ltmin=lt
+!     if(ldt.le.0.d0) ltmax=lt
+     ltn = lt+ldt
+     ltn = min(ltn,ltmax)
+     ltn = max(ltn,ltmin)
+     lt1=lt
+     lt=ltn
+     eps1=eps
+!     write(*,"(i4,1P12E19.10)") i,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0,d2,ldt
+     call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3)
+!     call findthis(lr,lt,y,eps,epst,d1,d2,d3)
+     if (abs(eps - eps0).lt.tol*abs(eps0)) then
+!        write(*,"(1P12E19.10)") tol,abs(eps-eps0)/eps0
+        exit
+     endif
+     !setup new d2
+
+     ! if we are closer than 10^-2  to the 
+     ! root (eps-eps0)=0, we are switching to 
+     ! the secant method, since the table is rather coarse and the
+     ! derivatives may be garbage.
+     if(abs(eps-eps0).lt.1.0d-3*abs(eps0)) then
+        d2 = (eps-eps1)/(lt-lt1)
+     endif
+!     if(i.ge.10) then
+!       write(*,*) "EOS: Did not converge in findtemp!"
+!       write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0"
+!     if(i.gt.5) then
+!       write(*,"(i4,1P10E22.14)") i,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0
+!     endif
+  enddo
+
+
+ if(i.ge.itmax) then
+    keyerrt=667
+    call bisection(lr,lt0,y,eps0,lt,alltables(:,:,:,2),keyerrt,1)
+    if(keyerrt.eq.667) then
+       if(lt.ge.log10(t_max_hack)) then
+          ! handling for too high temperatures
+          lt = log10(t_max_hack)
+          keyerrt=0
+          goto 12
+       else if(abs(lt-log10(t_max_hack))/log10(t_max_hack).lt.0.025d0) then
+          lt0 = min(lt,log10(t_max_hack))
+          keyerrt=0
+          goto 12
+       else
+          ! total failure
+          write(*,*) "EOS: Did not converge in findtemp!"
+          write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0"
+          write(*,"(i4,i4,1P10E19.10)") i,rl,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0
+          write(*,*) "Tried calling bisection... didn't help... :-/"
+          write(*,*) "Bisection error: ",keyerrt
+       endif
+    endif
+    
+    lt0=min(lt,log10(t_max_hack))
+    return
+ endif
+
+12 continue
+
+  lt0=min(lt,log10(t_max_hack))
+
+
+end subroutine findtemp
+
+subroutine findtemp_entropy(lr,lt0,y,sin,keyerrt,rfeps)
+
+! This routine finds the new temperature based
+! on rho, Y_e, entropy
+
+  use eosmodule
+
+  implicit none
+
+  real*8 lr,lt0,y,sin
+  real*8 s,lt,ldt
+  real*8 tol
+  real*8 d1,d2,d3
+  real*8 s0,s1,lt1
+
+  real*8 ltn,ltmax,ltmin
+  real*8 tinput,rfeps
+
+  integer :: rl = 0
+  integer itmax,i,keyerrt
+  integer ii,jj,kk
+
+  keyerrt=0
+
+  tol=rfeps ! need to find energy to less than 1 in 10^-10
+  itmax=20 ! use at most 20 iterations, then bomb
+
+  lt=lt0
+  lt1=lt 
+
+  s0=sin
+  s1=s0
+
+  ltmax=logtemp(ntemp)
+  ltmin=logtemp(1)
+
+  !preconditioning 1: do we already have the right temperature?
+  call findthis(lr,lt,y,s,alltables(:,:,:,3),d1,d2,d3)
+  if (abs(s-s0).lt.tol*abs(s0)) then
+     return
+  endif
+  lt1=lt
+  s1=s
+ 
+
+  do i=1,itmax
+     !d2 is the derivative ds/dlogtemp;
+     ldt = -(s - s0)/d2 
+     ltn = lt+ldt
+     ltn = min(ltn,ltmax)
+     ltn = max(ltn,ltmin)
+     lt1=lt
+     lt=ltn
+     s1=s
+     call findthis(lr,lt,y,s,alltables(:,:,:,3),d1,d2,d3)
+     if (abs(s - s0).lt.tol*abs(s0)) then
+       exit
+     endif
+     !setup new d2
+
+     ! if we are closer than 10^-2  to the 
+     ! root (eps-eps0)=0, we are switching to 
+     ! the secant method, since the table is rather coarse and the
+     ! derivatives may be garbage.
+     if(abs(s-s0).lt.1.0d-3*abs(s0)) then
+        d2 = (s-s1)/(lt-lt1)
+     endif
+  enddo
+
+
+ if(i.ge.itmax) then
+    keyerrt=667
+   call bisection(lr,lt0,y,s0,lt,alltables(:,:,:,3),keyerrt,2)
+    if(keyerrt.eq.667) then
+          write(*,*) "EOS: Did not converge in findtemp_entropy!"
+          write(*,*) "rl,logrho,logtemp0,ye,lt,s,s0,abs(s-s0)/s0"
+          write(*,"(i4,i4,1P10E19.10)") i,rl,lr,lt0,y,lt,s,s0,abs(s-s0)/s0
+          write(*,*) "Tried calling bisection... didn't help... :-/"
+          write(*,*) "Bisection error: ",keyerrt
+    endif
+    
+    lt0=lt
+    return
+ endif
+
+
+  lt0=lt
+
+
+end subroutine findtemp_entropy
diff --git a/src/nuc_eos/linterp_many.F90 b/src/nuc_eos/linterp_many.F90
new file mode 100644
index 0000000..05d9c7a
--- /dev/null
+++ b/src/nuc_eos/linterp_many.F90
@@ -0,0 +1,125 @@
+      SUBROUTINE intp3d_many ( x, y, z, f, kt, ft, nx, ny, nz, nvars, xt, yt, zt)
+!
+      implicit none
+!                                                          
+!---------------------------------------------------------------------
+!
+!     purpose: interpolation of a function of three variables in an
+!              equidistant(!!!) table.
+!
+!     method:  8-point Lagrange linear interpolation formula          
+!
+!     x        input vector of first  variable
+!     y        input vector of second variable
+!     z        input vector of third  variable
+!
+!     f        output vector of interpolated function values
+!
+!     kt       vector length of input and output vectors
+!
+!     ft       3d array of tabulated function values
+!     nx       x-dimension of table
+!     ny       y-dimension of table
+!     nz       z-dimension of table
+!     xt       vector of x-coordinates of table
+!     yt       vector of y-coordinates of table
+!     zt       vector of z-coordinates of table
+!
+!---------------------------------------------------------------------
+
+      integer kt,nx,ny,nz,iv,nvars
+      real*8 :: ft(nx,ny,nz,nvars)
+
+      real*8 x(kt),y(kt),z(kt),f(kt,nvars)
+      real*8 xt(nx),yt(ny),zt(nz)
+      real*8 d1,d2,d3
+!
+!
+      integer,parameter :: ktx = 1
+      real*8  fh(ktx,8,nvars), delx(ktx), dely(ktx), delz(ktx), &
+           a1(ktx,nvars), a2(ktx,nvars), a3(ktx,nvars), a4(ktx,nvars), &
+           a5(ktx,nvars), a6(ktx,nvars), a7(ktx,nvars), a8(ktx,nvars)
+
+      real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi
+      integer n,ix,iy,iz
+
+      IF (kt .GT. ktx)  STOP'***KTX**'
+!
+!
+!------  determine spacing parameters of (equidistant!!!) table
+!
+      dx    = (xt(nx) - xt(1)) / FLOAT(nx-1)
+      dy    = (yt(ny) - yt(1)) / FLOAT(ny-1)
+      dz    = (zt(nz) - zt(1)) / FLOAT(nz-1)
+!
+      dxi   = 1. / dx
+      dyi   = 1. / dy
+      dzi   = 1. / dz
+!
+      dxyi  = dxi * dyi
+      dxzi  = dxi * dzi
+      dyzi  = dyi * dzi
+!
+      dxyzi = dxi * dyi * dzi
+!
+!
+!------- loop over all points to be interpolated
+!
+      dO  n = 1, kt                                            
+!
+!------- determine location in (equidistant!!!) table 
+!                                                                  
+         ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi )
+         iy = 2 + INT( (y(n) - yt(1) - 1.e-10) * dyi )
+         iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi )
+!                                                     
+         ix = MAX( 2, MIN( ix, nx ) )
+         iy = MAX( 2, MIN( iy, ny ) )
+         iz = MAX( 2, MIN( iz, nz ) )
+!
+!         write(*,*) iy-1,iy,iy+1
+!
+!------- set-up auxiliary arrays for Lagrange interpolation
+!                                                                 
+         delx(n) = xt(ix) - x(n)
+         dely(n) = yt(iy) - y(n)
+         delz(n) = zt(iz) - z(n)
+!      
+         do iv = 1, nvars
+            fh(n,1,iv) = ft(ix  , iy  , iz, iv  )                             
+            fh(n,2,iv) = ft(ix-1, iy  , iz, iv  )                             
+            fh(n,3,iv) = ft(ix  , iy-1, iz, iv  )                             
+            fh(n,4,iv) = ft(ix  , iy  , iz-1, iv)                             
+            fh(n,5,iv) = ft(ix-1, iy-1, iz, iv  )                             
+            fh(n,6,iv) = ft(ix-1, iy  , iz-1, iv)                             
+            fh(n,7,iv) = ft(ix  , iy-1, iz-1, iv)                             
+            fh(n,8,iv) = ft(ix-1, iy-1, iz-1, iv)                             
+!              
+!------ set up coefficients of the interpolation polynomial and 
+!       evaluate function values 
+            !                                                    
+            a1(n,iv) = fh(n,1,iv)                             
+            a2(n,iv) = dxi   * ( fh(n,2,iv) - fh(n,1,iv) )       
+            a3(n,iv) = dyi   * ( fh(n,3,iv) - fh(n,1,iv) )       
+            a4(n,iv) = dzi   * ( fh(n,4,iv) - fh(n,1,iv) )       
+            a5(n,iv) = dxyi  * ( fh(n,5,iv) - fh(n,2,iv) - fh(n,3,iv) + fh(n,1,iv) )
+            a6(n,iv) = dxzi  * ( fh(n,6,iv) - fh(n,2,iv) - fh(n,4,iv) + fh(n,1,iv) )
+            a7(n,iv) = dyzi  * ( fh(n,7,iv) - fh(n,3,iv) - fh(n,4,iv) + fh(n,1,iv) )
+            a8(n,iv) = dxyzi * ( fh(n,8,iv) - fh(n,1,iv) + fh(n,2,iv) + fh(n,3,iv) + &
+                 fh(n,4,iv) - fh(n,5,iv) - fh(n,6,iv) - fh(n,7,iv) )
+!
+            f(n,iv)  = a1(n,iv) +  a2(n,iv) * delx(n)                         &
+                 +  a3(n,iv) * dely(n)                         &
+                 +  a4(n,iv) * delz(n)                         &
+                 +  a5(n,iv) * delx(n) * dely(n)               &
+                 +  a6(n,iv) * delx(n) * delz(n)               &
+                 +  a7(n,iv) * dely(n) * delz(n)               &
+                 +  a8(n,iv) * delx(n) * dely(n) * delz(n)     
+!
+         enddo
+
+      enddo
+!                                                                    
+      
+    end SUBROUTINE intp3d_many
+ 
diff --git a/src/nuc_eos/make.code.defn b/src/nuc_eos/make.code.defn
new file mode 100644
index 0000000..9be5ee3
--- /dev/null
+++ b/src/nuc_eos/make.code.defn
@@ -0,0 +1,7 @@
+SRCS = eosmodule.F90 nuc_eos.F90 bisection.F90 \
+       findtemp.F90 linterp_many.F90 readtable.F90
+
+SUBDIRS = 
+
+F90FLAGS+=-I$(HDF5_INC_DIRS) -I$(HDF5_LIB_DIRS)
+
diff --git a/src/nuc_eos/make.code.deps b/src/nuc_eos/make.code.deps
new file mode 100644
index 0000000..2b2549d
--- /dev/null
+++ b/src/nuc_eos/make.code.deps
@@ -0,0 +1,6 @@
+bisection.F90.o:       eosmodule.F90.o
+findtemp.F90.o:        eosmodule.F90.o
+linterp_many.F90.o:    eosmodule.F90.o
+nuc_eos.F90.o:         eosmodule.F90.o
+readtable.F90.o:       eosmodule.F90.o
+
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90
new file mode 100644
index 0000000..ba16882
--- /dev/null
+++ b/src/nuc_eos/nuc_eos.F90
@@ -0,0 +1,322 @@
+! #########################################################
+!
+! Copyright C. D. Ott and Evan O'Connor, July 2009
+!
+! UNITS: density              g/cm^3
+!        temperature          MeV
+!        ye                   number fraction per baryon
+!        energy               erg/g
+!        pressure             dyn/cm^2
+!        chemical potentials  MeV
+!        entropy              k_B / baryon
+!        cs2                  cm^2/s^2 (not relativistic)
+!
+! keyerr --> error output; should be 0
+! rfeps --> root finding relative accuracy, set around 1.0d-10
+! keytemp: 0 -> coming in with eps
+!          1 -> coming in with temperature
+!          2 -> coming in with entropy
+!
+subroutine nuc_eos_full(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
+     xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp,xabar,xzbar,xmu_e,xmu_n,xmu_p, &
+     xmuhat,keytemp,keyerr,rfeps)
+
+  use eosmodule
+  implicit none
+
+  real*8, intent(in)    :: xrho,xye
+  real*8, intent(inout) :: xtemp,xenr,xent
+  real*8, intent(out)   :: xprs,xcs2,xdedt
+  real*8, intent(out)   :: xdpderho,xdpdrhoe,xxa,xxh,xxn,xxp
+  real*8, intent(out)   :: xabar,xzbar,xmu_e,xmu_n,xmu_p,xmuhat
+  real*8, intent(in)    :: rfeps
+  integer, intent(in)   :: keytemp
+  integer, intent(out)  :: keyerr
+
+
+  ! local variables
+  real*8 :: lr,lt,y,xx,xeps,leps,xs
+  real*8 :: d1,d2,d3
+  real*8 :: ff(nvars)
+  integer :: keyerrt = 0
+
+  if(xrho.gt.eos_rhomax) then
+     stop "nuc_eos: rho > rhomax"
+  endif
+
+  if(xrho.lt.eos_rhomin) then
+     stop "nuc_eos: rho < rhomin"
+  endif
+
+  if(xye.gt.eos_yemax) then
+     stop "nuc_eos: ye > yemax"
+  endif
+
+  if(xye.lt.eos_yemin) then
+     stop "nuc_eos: ye < yemin"
+  endif
+
+  if(keytemp.eq.1) then
+     if(xtemp.gt.eos_tempmax) then
+        stop "nuc_eos: temp > tempmax"
+     endif
+     
+     if(xtemp.lt.eos_tempmin) then
+        stop "nuc_eos: temp < tempmin"
+     endif
+  endif
+
+  lr = log10(xrho)
+  lt = log10(xtemp)
+  y = xye
+  xeps = xenr + energy_shift
+  leps = log10(max(xeps,1.0d0))
+
+  keyerr = 0
+
+  if(keytemp.eq.0) then
+     !need to find temperature based on xeps
+     call findtemp(lr,lt,y,leps,keyerrt,rfeps)
+     if(keyerrt.ne.0) then
+        stop "Did not find temperature"
+     endif
+     xtemp = 10.0d0**lt
+
+  elseif(keytemp.eq.2) then
+     !need to find temperature based on xent
+     xs = xent
+     call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps)
+     xtemp = 10.0d0**lt
+  endif
+
+  ! have temperature, proceed:
+  call findall(lr,lt,y,ff)
+
+  xprs = 10.0d0**ff(1)
+
+  xeps = 10.0d0**ff(2) - energy_shift
+  if(keytemp.eq.1.or.keytemp.eq.2) then
+     xenr = xeps
+  endif
+
+  if(keytemp.eq.1.or.keytemp.eq.0) then
+     xent = ff(3)
+  endif
+
+  xcs2 = ff(5)
+
+! derivatives
+  xdedt = ff(6)
+
+  xdpdrhoe = ff(7)
+
+  xdpderho = ff(8)
+
+! chemical potentials
+  xmuhat = ff(9)
+
+  xmu_e = ff(10)
+
+  xmu_p = ff(11)
+
+  xmu_n = ff(12)
+
+! compositions
+  xxa = ff(13)
+
+  xxh = ff(14)
+
+  xxn = ff(15)
+
+  xxp = ff(16)
+
+  xabar = ff(17)
+
+  xzbar = ff(18)
+
+
+end subroutine nuc_eos_full
+
+subroutine nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+
+  implicit none
+  real*8 xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe
+  real*8,parameter :: idealK1 =  1.2435d15 * (0.5d0**(4.d0/3.d0))
+  real*8,parameter :: idealgamma = 1.41d0
+  integer keytemp
+
+  if(keytemp.eq.1) then
+!	energy wanted
+     xprs=idealK1*xrho**(idealgamma)
+     xenr=xprs/xrho/(idealgamma-1.d0)
+     xcs2=idealgamma*xprs/xrho
+  endif
+
+  xprs = (idealgamma - 1.0d0) * xrho * xenr
+
+  xdpderho = (idealgamma - 1.0d0 ) * xrho
+  xdpdrhoe = (idealgamma - 1.0d0 ) * xenr
+
+  xcs2= xdpdrhoe+xdpderho*xprs/xrho**2
+
+end subroutine nuc_low_eos
+
+subroutine nuc_eos_short(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
+     xdpderho,xdpdrhoe,xmunu,keytemp,keyerr,rfeps)
+
+  use eosmodule
+  implicit none
+
+  real*8, intent(in)    :: xrho,xye
+  real*8, intent(inout) :: xtemp,xenr,xent
+  real*8, intent(in)    :: rfeps
+  real*8, intent(out)   :: xprs,xmunu,xcs2,xdedt
+  real*8, intent(out)   :: xdpderho,xdpdrhoe
+  integer, intent(in)   :: keytemp
+  integer, intent(out)  :: keyerr
+
+  ! local variables
+  real*8 :: lr,lt,y,xx,xeps,leps,xs
+  real*8 :: d1,d2,d3,ff(8)
+  integer :: keyerrt = 0
+
+  if(xrho.gt.eos_rhomax) then
+     stop "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
+     return
+  endif
+
+  if(xye.gt.eos_yemax) then
+     stop "nuc_eos: ye > yemax"
+  endif
+
+  if(xye.lt.eos_yemin) then
+     stop "nuc_eos: ye < yemin"
+  endif
+
+  if(keytemp.eq.1) then
+     if(xtemp.gt.eos_tempmax) then
+        stop "nuc_eos: temp > tempmax"
+     endif
+     
+     if(xtemp.lt.eos_tempmin) then
+        call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp)
+        xent = 4.0d0
+        return
+     endif
+  endif
+
+  lr = log10(xrho)
+  lt = log10(xtemp)
+  y = xye
+  xeps = xenr + energy_shift
+  leps = log10(max(xeps,1.0d0))
+
+  keyerr = 0
+
+  if(keytemp.eq.0) then
+     !need to find temperature based on xeps
+     call findtemp(lr,lt,y,leps,keyerrt,rfeps)
+     if(keyerrt.ne.0) then
+        keyerr = keyerrt
+        return
+     endif
+     xtemp = 10.0d0**lt
+
+  elseif(keytemp.eq.2) then
+     !need to find temperature based on xent
+     xs = xent
+     call findtemp_entropy(lr,lt,y,xs,keyerrt,rfeps)
+     if(keyerrt.ne.0) then
+        keyerr = keyerrt
+        return
+     endif
+     xtemp = 10.0d0**lt
+  endif
+
+  ! have temperature, proceed:
+  call findall_short(lr,lt,y,ff)
+  xprs = 10.0d0**ff(1)
+
+  xeps = 10.0d0**ff(2) - energy_shift
+  if((keytemp.eq.1).or.(keytemp.eq.2)) then
+     xenr = xeps
+  endif
+
+  if(keytemp.eq.1.or.keytemp.eq.0) then
+     xent = ff(3)
+  endif
+
+  xmunu = ff(4)
+
+  xcs2 = ff(5)
+
+  xdedt = ff(6)
+
+  xdpdrhoe = ff(7)
+
+  xdpderho = ff(8)
+
+end subroutine nuc_eos_short
+
+subroutine findthis(lr,lt,y,value,array,d1,d2,d3)
+
+  use eosmodule
+
+  implicit none
+
+  integer rip,rim
+  integer tip,tim
+  integer yip,yim
+
+  real*8 lr,lt,y,value,d1,d2,d3
+  real*8 array(*)
+
+! Ewald's interpolator           
+  call intp3d(lr,lt,y,value,1,array,nrho,ntemp,nye,logrho,logtemp,ye,d1,d2,d3)
+
+
+end subroutine findthis
+
+
+subroutine findall(lr,lt,y,ff)
+
+  use eosmodule
+  implicit none
+
+  real*8 ff(nvars)
+  real*8 ffx(nvars,1)
+  real*8 lr,lt,y
+  integer i
+
+! Ewald's interpolator           
+  call intp3d_many(lr,lt,y,ffx,1,alltables,&
+       nrho,ntemp,nye,nvars,logrho,logtemp,ye)
+  ff(:) = ffx(:,1)
+
+
+end subroutine findall
+
+
+subroutine findall_short(lr,lt,y,ff)
+
+  use eosmodule
+  implicit none
+
+  real*8 ffx(8,1)
+  real*8 ff(8)
+  real*8 lr,lt,y
+  integer i
+  integer :: nvarsx = 8
+
+
+! Ewald's interpolator           
+  call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,1:8), &
+       nrho,ntemp,nye,nvarsx,logrho,logtemp,ye)
+  ff(:) = ffx(:,1)
+
+end subroutine findall_short
diff --git a/src/nuc_eos/readtable.F90 b/src/nuc_eos/readtable.F90
new file mode 100644
index 0000000..41f8159
--- /dev/null
+++ b/src/nuc_eos/readtable.F90
@@ -0,0 +1,234 @@
+subroutine readtable(eos_filename)
+! This routine reads the table and initializes
+! all variables in the module. 
+
+  use eosmodule
+  use hdf5 
+
+  implicit none
+
+  character(*) eos_filename
+
+  character(len=100) message
+
+! HDF5 vars
+  integer(HID_T) file_id,dset_id,dspace_id
+  integer(HSIZE_T) dims1(1), dims3(3)
+  integer error,rank,accerr
+  integer i,j,k
+
+  real*8 amu_cgs_andi
+  real*8 buffer1,buffer2,buffer3,buffer4
+  accerr=0
+
+  write(*,*) "Reading Ott EOS Table"
+
+  call h5open_f(error)
+
+  call h5fopen_f (trim(adjustl(eos_filename)), H5F_ACC_RDONLY_F, file_id, error)
+
+  write(6,*) trim(adjustl(eos_filename))
+
+! read scalars
+  dims1(1)=1
+  call h5dopen_f(file_id, "pointsrho", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_INTEGER, nrho, dims1, error)
+  call h5dclose_f(dset_id,error)
+
+  if(error.ne.0) then
+     stop "Could not read EOS table file"
+  endif
+
+  dims1(1)=1
+  call h5dopen_f(file_id, "pointstemp", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_INTEGER, ntemp, dims1, error)
+  call h5dclose_f(dset_id,error)
+
+  if(error.ne.0) then
+     stop "Could not read EOS table file"
+  endif
+
+  dims1(1)=1
+  call h5dopen_f(file_id, "pointsye", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_INTEGER, nye, dims1, error)
+  call h5dclose_f(dset_id,error)
+
+  if(error.ne.0) then
+     stop "Could not read EOS table file"
+  endif
+
+  write(message,"(a25,1P3i5)") "We have nrho ntemp nye: ", nrho,ntemp,nye
+  write(*,*) message
+
+  allocate(alltables(nrho,ntemp,nye,nvars))
+
+  ! index variable mapping:
+  !  1 -> logpress
+  !  2 -> logenergy
+  !  3 -> entropy
+  !  4 -> munu
+  !  5 -> cs2
+  !  6 -> dedT
+  !  7 -> dpdrhoe
+  !  8 -> dpderho
+  !  9 -> muhat
+  ! 10 -> mu_e
+  ! 11 -> mu_p
+  ! 12 -> mu_n
+  ! 13 -> xa
+  ! 14 -> xh
+  ! 15 -> xn
+  ! 16 -> xp
+  ! 17 -> abar
+  ! 18 -> zbar
+
+
+  dims3(1)=nrho
+  dims3(2)=ntemp
+  dims3(3)=nye
+  call h5dopen_f(file_id, "logpress", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,1), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+  call h5dopen_f(file_id, "logenergy", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,2), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+  call h5dopen_f(file_id, "entropy", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,3), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+  call h5dopen_f(file_id, "munu", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,4), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+  call h5dopen_f(file_id, "cs2", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,5), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+  call h5dopen_f(file_id, "dedt", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,6), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+  call h5dopen_f(file_id, "dpdrhoe", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,7), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+  call h5dopen_f(file_id, "dpderho", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,8), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+! chemical potentials
+  call h5dopen_f(file_id, "muhat", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,9), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  call h5dopen_f(file_id, "mu_e", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,10), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  call h5dopen_f(file_id, "mu_p", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,11), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  call h5dopen_f(file_id, "mu_n", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,12), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+! compositions
+  call h5dopen_f(file_id, "Xa", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,13), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  call h5dopen_f(file_id, "Xh", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,14), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  call h5dopen_f(file_id, "Xn", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,15), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  call h5dopen_f(file_id, "Xp", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,16), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+
+! average nucleus
+  call h5dopen_f(file_id, "Abar", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,17), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  call h5dopen_f(file_id, "Zbar", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,18), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+! Gamma
+  call h5dopen_f(file_id, "gamma", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, alltables(:,:,:,19), dims3, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  allocate(logrho(nrho))
+  dims1(1)=nrho
+  call h5dopen_f(file_id, "logrho", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, logrho, dims1, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  allocate(logtemp(ntemp))
+  dims1(1)=ntemp
+  call h5dopen_f(file_id, "logtemp", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, logtemp, dims1, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  allocate(ye(nye))
+  dims1(1)=nye
+  call h5dopen_f(file_id, "ye", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, ye, dims1, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+
+  call h5dopen_f(file_id, "energy_shift", dset_id, error)
+  call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, energy_shift, dims1, error)
+  call h5dclose_f(dset_id,error)
+  accerr=accerr+error
+
+  if(accerr.ne.0) then
+    stop "Problem reading EOS table file"
+  endif
+
+
+  call h5fclose_f (file_id,error)
+
+  call h5close_f (error)
+
+  ! set min-max values:
+
+  eos_rhomin = 10.0d0**logrho(1)
+  eos_rhomax = 10.0d0**logrho(nrho)
+
+  eos_yemin = ye(1)
+  eos_yemax = ye(nye)
+
+  eos_tempmin = 10.0d0**logtemp(1)
+  eos_tempmax = 10.0d0**logtemp(ntemp)
+
+  write(6,*) "Done reading eos tables"
+
+
+end subroutine readtable
+
+