anchor.diff: -! a simple gamma-law equation of state
  -! p = rho e (gamma - 1)
  +! a simple gamma-law equation of state: 
  +!    p = rho e (gamma - 1)
  +!
  +! to define the temperature and entropy, we will use the ideal gas law:
  +!    p = rho k T / (mu m_amu)
  +!    
  +! where mu is the mean molecular weight and m_amu is the atomic mass
  +! unit
  +
   module eos_module
   
     use datatypes_module
     use params_module, only: gamma
     implicit none
   
  -  integer, parameter :: eos_input_p = 1
  -  integer, parameter :: eos_input_e = 2
  +  integer, parameter :: eos_input_p  = 1
  +  integer, parameter :: eos_input_e  = 2
     integer, parameter :: eos_input_pe = 3
  -
  +  integer, parameter :: eos_input_ps = 4
     private
   
  -  public :: eos_input_p, eos_input_e, eos_input_pe
  +  public :: eos_input_p, eos_input_e, eos_input_pe, eos_input_ps
     public :: eos
   
   contains
   
  -  subroutine eos(input, p, e, rho)
  +  subroutine eos(input, p, e, rho, T, s)
       
       integer,          intent(in   ) :: input
       real (kind=dp_t), intent(inout) :: p
       real (kind=dp_t), intent(inout) :: e
       real (kind=dp_t), intent(inout) :: rho
  +    real (kind=dp_t), intent(inout), optional :: s
  +    real (kind=dp_t), intent(inout), optional :: T
  +
  +    ! physical constants
  +    real (kind=dp_t), parameter :: mu = 1.0_dp_t
  +    real (kind=dp_t), parameter :: m_amu = 1.66053886e-24_dp_t   ! g
  +    real (kind=dp_t), parameter :: k = 1.3806503e-16_dp_t        ! erg / K
  +
  +    real (kind=dp_t), parameter :: c_v = 1.5_dp_t*k/(mu*m_amu)
   
       if (input == eos_input_p) then
   
  
positive.code: 
         rho = p/e/(gamma - 1.0_dp_t)
  
      else
         print *, 'ERROR: invalid input in EOS'
         stop
      endif
  
    end subroutine eos
  
  end module eos_module
  
negative.code: 
         e = p / (rho * (gamma - 1.0_dp_t))
  
         if (present(T)) T = mu*m_amu*e/(1.5_dp_t*k)
         if (present(s)) s = c_v*log(p/rho**gamma) 
  
      else if (input == eos_input_e) then
  
         if (e < 0.0_dp_t) then
  