anchor.diff:          rpar(iTeff) = Teff
           rpar(itau) = tau
           ipar(ihandle) = eos_handle
  -
  +        lrpar = number_photosphere_rpar
  +        lipar = number_photosphere_ipar
  +        
           ! 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
  +        sigma_Th = 8.0_dp*onethird*pi*(electroncharge**2/Melectron/clight2)**2
  +        kappa_Th = sigma_Th* Z/A
  +        fallback_Pphoto = 2.0_dp*onethird*arad*Teff**4   
  +        Pgas = twothird*grav/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   
  +        rho_guess = Pgas*amu*A/(Z+1.0)/(boltzmann*Teff)
           ! get brackets for root find
           drho = 0.1_dp
           maximum_iterations = 10
  
positive.code:     
      subroutine find_photospheric_pressure(Teff,grav,tau,Pphoto,eos_handle,ierr)
          use constants_def
          use nucchem_def, only : nuclide_not_found
      	use nucchem_lib, only : get_nuclide_index
      	use num_lib
          integer, parameter :: default_maximum_iterations_photosphere = 20
  
negative.code:         real(dp) :: root_ph,rho,dfdrho
          integer, pointer :: ipar(:) => null() ! (lipar)
          real(dp), pointer :: rpar(:) => null()  ! (lrpar)
          integer :: i, lrpar, lipar
          real(dp) :: sigma_Th
          real(dp) :: fallback_Pphoto
          real(dp) :: rho_guess, kappa_Th, A, Z, Pgas
          real(dp) :: rho1, rho3, drho, y1, y3
  