anchor.diff:          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)
  -
  +
  +        ! set size of data structure        
  +        lrpar = number_base_rpar + number_species
  +        lipar = number_base_ipar + number_species
  +        allocate(ipar(lipar), rpar(lrpar))
  +
  +        ! set photosphere to be pure He and compute moments
  +        Y = [1.0_dp, 0.0_dp]/nuclib% A(chem_ids)
           call compute_composition_moments(number_species,chem_ids,Y,ionic,Xsum, &
           &   ncharged, charged_ids, Yion, exclude_neutrons = .TRUE.)
   
  +        ! sanity check on composition
           if (Xsum - 1.0_dp > 2.0*epsilon(1.0_dp)) then
               ierr = bad_composition
               return
  
positive.code: module bc09
      
      ! data structure for the parameter arrays
      ! moments of the composition
  
negative.code:         maximum_iterations = default_maximum_iterations_photosphere
          eps_lnrho = default_tolerance_photosphere_lnrho
          eps_ph = default_tolerance_photosphere_condition
          fallback_Pphoto = 2.0_dp*onethird*arad*Teff**4
  
          ! scale tolerance to a thomson scaterring atmosphere
          sigma_Th = 8.0_dp*onethird*pi*(electroncharge**2/Melectron/clight2)**2
  