anchor.diff:      use dStar_eos_lib
       use conductivity_lib
       use PPP_electron
  -    use num_lib, only: binary_search
       
       implicit none
  -    
  +    integer, parameter :: NX=11,Nrho=13,NT=5
  +    real(dp), parameter :: Xmin=0.0_dp, Xmax = 1.0_dp
  +    real(dp), parameter :: lgrho_min=5.0_dp, lgrho_max = 9.0_dp
  +    real(dp), parameter :: lgT_min=7.0_dp, lgT_max = 9.0_dp
       integer :: eos_handle,cond_handle,ierr,i,j,k
       integer :: Z(2), N(2)
       real(dp) :: A(2), X(2)
  
positive.code:         &   (trim(name(j)),X(j),j=1,2)
          write(output_unit,'(2(a,f6.2,tr4))') '<Z> = ',ionic%Z,'<A> = ',ionic% A
          write (output_unit, &
          &   '(2a6,6a11,a7/,2("======"),6("==========="),"=======")') &
              & 'lg(r)','lg(T)','Gamma','eta_e', &
              & 'K_ee','K_ei','K_tot','K_table','diff'
              
          temperature: do j = 1, 5
              lgT = real(j-1,dp)*0.5_dp + 7.0_dp
              T = 10.0_dp**lgT
              
              density: do i = 1, 13
                  lgrho = real(i-1,dp)/3.0_dp + 5.0_dp
                  rho = 10.0_dp**lgrho
                  chi = use_default_nuclear_size
                  call eval_crust_eos(eos_handle,rho,T,ionic, &
                      & ncharged, charged_ids, Yion, Tcs,  &
  
negative.code:     call check_okay% assert(ierr==0)
  
      diff = 0.0_dp
      call conductivity_set_controls(cond_handle,which_ee_scattering=icond_pcy)    
      Tcs = 1.0e9_dp
      Z = [ 26, 2 ]
      N = [ 30, 2 ]
      A = real(Z+N,dp)
      chem_ids = [ (get_nuclide_index_from_ZN(Z(k),N(k)),k=1,2) ]
      name = [ (nuclib% name(chem_ids(j)), j=1,2) ]
      name(1:2)(1:1) = [ (StrUpCase(name(j)(1:1)), j=1,2) ]
      Z1 = real(Z(1),dp); Z2 = real(Z(2),dp); A1 = real(A(1),dp); A2 = real(A(2),dp)
  
      X1 = linspace(NX,Xmin,Xmax)
      lgrho = linspace(Nrho,lgrho_min,lgrho_max)
      lgT = linspace(NT,lgT_min,lgT_max)
      
      composition: do k = 1, NX
          X(1) = X1(k)
          X(2) = 1.0_dp -X(1)
          Y = X/A
          call compute_composition_moments(2,chem_ids,Y,ionic,Xsum, &
  