anchor.diff:    subroutine variables(T, Sigma, Omega, f, optical_depth)
       implicit none
   
  -    real(kind = x_precision),intent(in)          :: T,Sigma,Omega
  -    integer, intent(in)                          :: optical_depth
  -    real(kind = x_precision),intent(out)         :: f
  -    real(kind = x_precision)                     :: coeff_a=0.,coeff_b=0.,coeff_c=0.
  -
  -    real(kind = x_precision)                     :: H
  -    real(kind = x_precision)                     :: rho
  -    real(kind = x_precision)                     :: cs
  -    real(kind = x_precision)                     :: nu
  -    real(kind = x_precision)                     :: Q_plus
  -    real(kind = x_precision)                     :: K_ff
  -    real(kind = x_precision)                     :: K_e
  -    real(kind = x_precision)                     :: E_ff
  -    real(kind = x_precision)                     :: tau_eff
  -    real(kind = x_precision)                     :: Fz
  -    real(kind = x_precision)                     :: Q_minus
  +    real(kind = x_precision),intent(in)  :: T,Sigma,Omega
  +    integer, intent(in)                  :: optical_depth
  +    real(kind = x_precision),intent(out) :: f
  +    real(kind = x_precision)             :: coeff_a=0.,coeff_b=0.,coeff_c=0.
  +
  +    real(kind = x_precision)             :: H
  +    real(kind = x_precision)             :: rho
  +    real(kind = x_precision)             :: cs
  +    real(kind = x_precision)             :: nu
  +    real(kind = x_precision)             :: Q_plus
  +    real(kind = x_precision)             :: K_ff
  +    real(kind = x_precision)             :: K_e
  +    real(kind = x_precision)             :: E_ff
  +    real(kind = x_precision)             :: tau_eff
  +    real(kind = x_precision)             :: Fz
  +    real(kind = x_precision)             :: Q_minus
       !------------------------------------------------------------------------
   
       coeff_a = (Omega * state_0%Omega_0)**2 * Sigma * state_0%S_0 / 2._x_precision
  
positive.code:   real(kind=x_precision) function dichotomy(Smin, Smax, eps, T, omega, optical_depth)
      implicit none
  
      real(kind=x_precision),intent(inout)                     :: Smin,Smax
      real(kind=x_precision),intent(in)                        :: eps
      real(kind=x_precision),intent(in)                        :: T
      real(kind=x_precision),intent(in)                        :: omega
      integer,intent(in)                                       :: optical_depth
  
      real(kind=x_precision)                                   :: f_min
      real(kind=x_precision)                                   :: f_max
      real(kind=x_precision)                                   :: f_center
      real(kind=x_precision)                                   :: S_center = 0._x_precision
      integer                                                  :: j = 0
  
      !-------------------------------------------------------------------------
      ! N-> Number of iterations for the dichotomy
  
negative.code:          index_fcp, index_scp, sigma_c_thin, temp_c_thin, nb_it)
      implicit none
  
      integer                                  ,intent(in)  :: nb_it
      real(kind = x_precision),dimension(nb_it),intent(in)  :: sigma_real_thick
      real(kind = x_precision),dimension(nb_it),intent(in)  :: temp_real_thick
      real(kind = x_precision),dimension(nb_it),intent(in)  :: sigma_real_thin
  
      integer                                  ,intent(in)  :: index_fcp
      integer                                  ,intent(out) :: index_scp
  
      real(kind = x_precision)                 ,intent(out) :: sigma_c_thin
      real(kind = x_precision)                 ,intent(out) :: temp_c_thin
      integer                                               :: i
      i = max(1, index_fcp)
      
      ! The change occurs when the difference of sigma changes sign
  