anchor.diff:       ! Use the ambient gas composition (we don't need to worry about saving
        ! our result in eint_ambient since it will be reset in the next call).
   
  -     call eos_given_RTX(eint_ambient,small_pres,small_dens, &
  -                        small_temp,xn_ambient)
  +     eos_state % xn  = xn_ambient
  +     eos_state % rho = small_dens
  +     eos_state % T   = small_temp
  + 
  +     call eos(eos_input_rt, eos_state, .false.)
  +
  +     small_pres = eos_state % p
   
        ! Get the rest of the ambient thermodynamic state
  -     call eos_given_RTX(eint_ambient,pres_ambient,dens_ambient, &
  -                        temp_ambient,xn_ambient)
  +
  +     eos_state % rho = dens_ambient
  +     eos_state % T   = temp_ambient
  +
  +     call eos(eos_input_rt, eos_state, .false.)
  +
  +     eint_ambient = eos_state % e
  +     pres_ambient = eos_state % p
   
        ! Compute the radius of the primary
        radius_P_initial = -1.0d0
  
positive.code:      double precision :: pres_zone, temp_zone
       double precision :: dist_P, dist_S
  
       integer :: i,j,k,ii,jj,kk,n
  
       !$OMP PARALLEL DO PRIVATE(i, j, k, xl, yl, zl, xx, yy, zz) &
  
negative.code:      double precision :: a
       double precision :: length
  
       type (eos_t) :: eos_state
  
       ! Temporary storage variables in case we need to switch the primary and secondary.
  
       double precision, allocatable :: temp_model_state(:,:), temp_model_r(:)
  