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 ! precision to which we need to find temp
  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
 
  do i=1,itmax
     !d2 is the derivative deps/dlogtemp;
     ldt = -(eps - eps0)/d2 
     ltn = lt+ldt
     ltn = min(ltn,ltmax)
     ltn = max(ltn,ltmin)
     lt1=lt
     lt=ltn
     eps1=eps
     call findthis(lr,lt,y,eps,alltables(:,:,:,2),d1,d2,d3)
     if (abs(eps - eps0).lt.tol*abs(eps0)) 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(eps-eps0).lt.1.0d-3*abs(eps0)) then
        d2 = (eps-eps1)/(lt-lt1)
     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
#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
    
    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 0
    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
#endif
    
    lt0=lt
    return
 endif


  lt0=lt


end subroutine findtemp_entropy


subroutine findtemp_low(lr,t0,y,epsin,keyerrt,rfeps)

  ! this routine is for finding fake temperatures
  ! outside of the table range; we use linear 
  ! extrapolation in temperature

  use eosmodule

  implicit none

  real*8 lr,t0,y,epsin
  real*8 eps,t,dt
  real*8 tol
  real*8 dlepsdt
  real*8 eps0,eps1,t1

  real*8 tn,tmin,tmax
  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=100 ! use at most 20 iterations, then bomb

  t=t0
  t1=t 

  eps0=epsin
  eps1=eps0
  
  tmax=10.0d0**logtemp(2)
  tmin=-20.0d0

  !preconditioning 1: do we already have the right temperature?
  call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,&
       ntemp,nye,logrho,logtemp,ye,dlepsdt)

  if (abs(eps-eps0).lt.tol*abs(eps0)) then
     return
  endif
  t1=t
  eps1=eps
 
  do i=1,itmax
     dt = -(eps - eps0)/dlepsdt 
     tn = t+dt
     if(tn >= tmax) then
        tn = tmax*0.9d0
     endif
     tn = min(tn,tmax)
     tn = max(tn,tmin)
     t1=t
     t=tn
     eps1=eps

     call intp3d_linearTlow(lr,t,y,eps,1,alltables(:,:,:,2),nrho,&
          ntemp,nye,logrho,logtemp,ye,dlepsdt)

     if (abs(eps - eps0).lt.tol*abs(eps0)) then
        exit
     endif

     ! 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
        dlepsdt = (eps-eps1)/(t-t1)
     endif
  enddo

  t0 = t

  if(i.ge.itmax) then
     keyerrt=667
     return
  endif


end subroutine findtemp_low