anchor.diff:         Yion=>rpar(number_base_rpar+1:number_base_rpar+ncharged)
          chi = use_default_nuclear_size
   
  +       call get_rho_from_PT()
  +       rpar(irho) = exp(lnrho)
          call eval_crust_eos(eos_handle,rho,Teff,ionic,ncharged,charged_ids,Yion, &
                  &   res,phase,chi)
   
  
positive.code:         ierr = 0
          
          ! composition is a He/Fe mix
          ! set size of data structure for rpar, ipar
          
          lrpar = number_base_rpar + number_species
          lipar = number_base_ipar + number_species
          allocate(ipar(lipar), rpar(lrpar))
  
          chem_ids = [ get_nuclide_index('he4'), get_nuclide_index('fe56') ]
          if (any(chem_ids == nuclide_not_found)) then    ! this is a fatal error
              ierr = nuclide_not_found
              return
          endif
          
          ! set photosphere to be pure He
          Y = [1.0_dp, 0.0_dp]
          Y = Y/nuclib% A(chem_ids)
  
          call compute_composition_moments(number_species,chem_ids,Y,ionic,Xsum, &
          &   ncharged, charged_ids, Yion, exclude_neutrons = .TRUE.)
  
          if (Xsum - 1.0_dp > 2.0*epsilon(1.0_dp)) then
              ierr = bad_composition
              return
  
negative.code:     integer, parameter :: iTeff = itau + 1
      integer, parameter :: itemp = iTeff + 1
      integer, parameter :: ipres = itemp + 1
      integer, parameter :: irho = ipres + 1
      integer, parameter :: iKph = irho + 1
      integer, parameter :: iChi_rho = iKph + 1
      integer, parameter :: iChi_T = iChi_rho + 1
      integer, parameter :: del_ad = iChi_T + 1
  