anchor.diff:      real(kind = x_precision)             :: Q_minus
       !------------------------------------------------------------------------
   
  +    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)
  
positive.code:     real(kind = x_precision)                               :: Smin
      real(kind = x_precision)                               :: Smax
  
      real(kind = x_precision)                               :: omega          ! Angular velocity
      real(kind = x_precision)                               :: f              ! Q+ - Q-
      real(kind = x_precision)                               :: tau_eff        ! Effective optical depth
      integer                                                :: optical_depth  ! Indicator for the optical thickness
  
negative.code:     else if( f_max * f_min < 0.) then
        iteration:do while (dabs(Smax - Smin) >= eps .and. j < max_it)
  
          call variables(k, T, Smin, f_min, optical_depth, tau_eff)
  
          call variables(k, T, Smax, f_max, optical_depth, tau_eff)
  
          call variables(k, T, S_center, f_center, optical_depth, tau_eff)
  
          if (f_min * f_center > 0.) then
  
  