2 files changed, 108 insertions, 2 deletions
diff --git a/src/nuc_eos/findtemp.F90 b/src/nuc_eos/findtemp.F90
index d208a5f..646994d 100644
--- a/src/nuc_eos/findtemp.F90
+++ b/src/nuc_eos/findtemp.F90
@@ -16,10 +16,10 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,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
+  tol=rfeps ! need to find energy to less than 1 in 10^10
   itmax=20 ! use at most 20 iterations, then bomb
 
   lt=lt0
@@ -95,12 +95,14 @@ subroutine findtemp(lr,lt0,y,epsin,keyerrt,rfeps)
           keyerrt=0
           goto 12
        else
+#if 0
           ! 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
     endif
     
diff --git a/src/nuc_eos/nuc_eos.F90 b/src/nuc_eos/nuc_eos.F90
index e84ac82..399d89b 100644
--- a/src/nuc_eos/nuc_eos.F90
+++ b/src/nuc_eos/nuc_eos.F90
@@ -348,6 +348,91 @@ subroutine nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,&
 
 end subroutine nuc_eos_press_eps
 
+subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,&
+     xdpderho,&
+     keytemp,keyerr,rfeps)
+
+  use eosmodule
+  implicit none
+
+  real*8, intent(in)    :: xrho,xye
+  real*8, intent(inout) :: xtemp,xenr
+  real*8, intent(in)    :: rfeps
+  real*8, intent(out)   :: xdpdrhoe,xdpderho
+  integer, intent(in)   :: keytemp
+  integer, intent(out)  :: keyerr
+
+  ! local variables
+  real*8 :: xcs2,xprs
+  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)
+     return
+  endif
+
+  if(xye.gt.eos_yemax) then
+     stop "nuc_eos: ye > yemax"
+  endif
+
+  if(xye.lt.eos_yemin) then
+     keyerr = 102
+     write(6,*) "ye: ", xye
+     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)
+        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
+     stop "eos_nuc_press does not support keytemp.eq.2"
+  endif
+
+  ! have temperature, proceed:
+  call findall_dpdr_dpde(lr,lt,y,ff)
+  xdpdrhoe = 10.0d0**ff(1)
+  xdpderho = 10.0d0**ff(2)
+
+
+  if((keytemp.eq.1).or.(keytemp.eq.2)) then
+     call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3)
+     xeps = 10.0d0**xeps - energy_shift
+     xenr = xeps
+  endif
+
+end subroutine nuc_eos_dpdr_dpde
+
 
 subroutine findthis(lr,lt,y,value,array,d1,d2,d3)
 
@@ -426,3 +511,22 @@ subroutine findall_press_eps(lr,lt,y,ff)
   ff(:) = ffx(:,1)
 
 end subroutine findall_press_eps
+
+subroutine findall_dpdr_dpde(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 = 2
+
+
+! Ewald's interpolator           
+  call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,7:8), &
+       nrho,ntemp,nye,nvarsx,logrho,logtemp,ye)
+  ff(:) = ffx(:,1)
+
+end subroutine findall_dpdr_dpde