anchor.diff:         rpar(ipres) = P
          rpar(itemp) = T
   
  -       call get_rho_from_PT()
          eta = res(i_Theta) !1.0/TpT
          Gamma = res(i_Gamma)
          call get_thermal_conductivity(rho,T,chi, &
  
positive.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
  
negative.code:        rpar(iKph) = kappa
  
         dlnT4dlnP(1) = 0.75_dp*kappa*P/grav/exp(lnT4(1))
         
      contains
          subroutine get_rho_from_PT(ierr)
              use num_lib, only: safe_root_without_brackets
              integer, intent(out) :: ierr
              lnrho = safe_root_without_brackets(eval_pressure,rho_guess,dlnrho,max_newt, &
              &   imax, epsx,epsy,lrpar,rpar,lipar,ipar,ierr)
              
          end subroutine get_rho_from_PT
  
      end subroutine deriv
  
      real(dp) function eval_pressure(lnrho, chi_rho, lrpar, rpar, lipar, ipar, ierr)
  