anchor.diff:      !------------------------------------------------------------------------
   
       Omega   = x_state%Omega(k)
  +
       coeff_a = (Omega * state_0%Omega_0)**2 * Sigma * state_0%S_0 / 2._x_precision
       coeff_b = (-1._x_precision/3._x_precision) * cst_rad * (T * state_0%T_0)**4 / state_0%H_0
       coeff_c = - params%RTM * T * Sigma * state_0%S_0 / (2._x_precision * state_0%H_0**2)
  -
       delta   = coeff_b**2 - 4._x_precision * coeff_a * coeff_c
       H       = -0.5_x_precision * (coeff_b + sign(sqrt(delta),coeff_b)) / coeff_a
   
  
positive.code:       temp_thin(k)   =  temp_c_thin * state_0%T_0
        sigma_thin(k)  =  sigma_c_thin * state_0%S_0
  
        call variables(k, temp_c_thick, sigma_c_thick, f, optical_depth, tau_eff)
        tau_thick(k)   = tau_eff
        call variables(k, temp_c_thin, sigma_c_thin, f, optical_depth, tau_eff)
  
negative.code:       dichotomy = 0
  
      else if( f_max * f_min < 0.) then
        do while (dabs(Smax - Smin) >= eps .and. j < max_it)
  
          call variables(k, T, Smin, f_min, optical_depth, tau_eff)
  
  