anchor.diff:  		real(dp) :: eta, g14
           integer ::  size_tab ! = 4*size(Tb)
           real(dp), dimension(:), allocatable :: tabTb9, tabTeff, tabTeff6_4
  +        integer :: i
           
           ! make a very dense table of Tb(Teff); then interpolate to get Teff(Tb)        
           size_tab = 4*size(Tb)
  
positive.code:         rpar(iTeff) = Teff
          rpar(itau) = tau
          ipar(ihandle) = eos_handle
  
          ! use initial guess with ideal gas pressure and thomson scattering
          A = nuclib% A(ipar(ichem_id))
          Z = nuclib% Z(ipar(ichem_id))
          kappa_Th = sigma_Th*nuclib% Z/A
          Pgas = twothird*gravity/kappa_th - onethird*arad*Teff**4
          if (Pgas < 0.0_dp) then
              ierr = negative_photosphere_gas_pressure
              Pphoto = fallback_Pphoto
              return
          end if
          rhoph_guess = Pgas*amu*A/(Z+1.0)/(boltzmann*Teff)
          fallback_Pphoto = 2.0_dp*onethird*arad*Teff**4   
          ! get brackets for root find
          drho = 0.1_dp
          maximum_iterations = 10
  
negative.code:         eps_rho = default_tolerance_photosphere_rho
          eps_ph = default_tolerance_photosphere_condition
  
  		root_ph = safe_root_with_initial_guess(photosphere,rho_guess,rho1,rho3,y1,y3, &
              &   maximum_iterations,eps_rho,eps_ph,lrpar,rpar,lipar,ipar,ierr)
  
          if (ierr /= 0) then
              write(*,*) 'unable to converge', rho_guess, rho1, rho3, y1, y3
              return
          end if
          Pphoto = rpar(iPph)
  