anchor.diff:         endif
   
          p = rho*e*(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_pe) then
   
  
positive.code: 
         e = p / (rho * (gamma - 1.0_dp_t))
  
      else if (input == eos_input_e) then
  
         if (e < 0.0_dp_t) then
  
negative.code: 
         rho = p/e/(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_ps) then
  
         ! entropy and pressure and inputs, return density and e 
         if (.not. present(s)) then
            print *, 'ERROR: entropy must be input'
            stop
         endif
         
         ! note: our definition of entropy neglects the constant since
         ! we really only care about changes in entropy.
         rho = (p*exp(-s/c_v))**(1.0_dp_t/gamma)
  
         e = p / (rho * (gamma - 1.0_dp_t))
  
         if (present(T)) T = mu*m_amu*e/(1.5_dp_t*k)
  
      else
         print *, 'ERROR: invalid input in EOS'
         stop
      endif
  
  
    end subroutine eos
  
  end module eos_module
  