anchor.diff:       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(:)
  
positive.code:      ! 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)
  
       ! Get the rest of the ambient thermodynamic state
       call eos_given_RTX(eint_ambient,pres_ambient,dens_ambient, &
                          temp_ambient,xn_ambient)
  
       ! Compute the radius of the primary
       radius_P_initial = -1.0d0
  
negative.code:               temp_zone = temp_zone/(nsub*nsub*nsub)
                state(i,j,k,UFS:UFS-1+nspec) = state(i,j,k,UFS:UFS-1+nspec)/(nsub*nsub*nsub)
  
                ! thermodynamics              
  
                eos_state % xn  = state(i,j,k,UFS:UFS-1+nspec)
                eos_state % p   = pres_zone
                eos_state % rho = state(i,j,k,URHO)
                eos_state % T   = temp_zone
  
                pt_index(1) = i
                pt_index(2) = j
                pt_index(3) = k
  
                if (interp_temp) then
                   state(i,j,k,UTEMP) = eos_state % T
                   call eos(eos_input_rt, eos_state, .false., pt_index = pt_index)
                else
                   call eos(eos_input_rp, eos_state, .false., pt_index = pt_index)
                   state(i,j,k,UTEMP) = eos_state % T
                endif
  
                state(i,j,k,UEINT) = eos_state % e
  
                state(i,j,k,UEDEN) = state(i,j,k,URHO) * state(i,j,k,UEINT)
                state(i,j,k,UEINT) = state(i,j,k,URHO) * state(i,j,k,UEINT)
  