anchor.diff:      integer                                                :: i     = 0
       integer                                                :: j     = 0
       real(kind = x_precision)                               :: temp  = 0.0d0
  -    real(kind = x_precision), parameter                    :: t_min = 1.2d-2
  -    real(kind = x_precision), parameter                    :: t_max = 1.0d0
  +    real(kind = x_precision), parameter                    :: t_min = 8.5d-3
  +    real(kind = x_precision), parameter                    :: t_max = 2.3d2
       integer,                  parameter                    :: nb_it = 10
       real(kind = x_precision)                               :: sigma = 0.0d0
       real(kind = x_precision)                               :: Smin  = 0.0d0
       real(kind = x_precision)                               :: Smax  = 0.0d0
  -    real(kind = x_precision)                               :: eps   = 1.0d-4
  +    real(kind = x_precision)                               :: eps   = 1.0d-9
       real(kind = x_precision)                               :: omega = 0.0d0
       real(kind = x_precision)                               :: r     = 5.0d0
       real(kind = x_precision)                               :: H     = 0.0d0    
       real(kind = x_precision)                               :: rho   = 0.0d0       
       real(kind = x_precision)                               :: cs    = 0.0d0
       real(kind = x_precision)                               :: nu    = 0.0d0  
  -    real(kind = x_precision)                               :: Q_plus = 0.0d0        
  +    real(kind = x_precision)                               :: Q_plus = 0.0d0    
       real(kind = x_precision)                               :: Q_minus= 0.0d0
       real(kind = x_precision)                               :: K_ff   = 0.0d0    
  -    real(kind = x_precision)                               :: K_e    = 0.0d0      
  +    real(kind = x_precision)                               :: K_e    = 0.0d0    
       real(kind = x_precision)                               :: tau_eff= 0.0d0
       real(kind = x_precision)                               :: P_rad  = 0.0d0    
       real(kind = x_precision)                               :: P_gaz  = 0.0d0
  
positive.code:   !K_ff, K_e, tau_eff, P_rad, P_gaz,E_ff,Fz for T, Sigma and Omega given
    !------------------------------------------------------------------------
    subroutine variables(T, Sigma, Omega, H, rho, cs, nu, Q_plus, Q_minus,&
         K_ff, K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f,Sigma_0)
      implicit none
   
      type(parameters)                                         :: param
      real(kind = x_precision),intent(in)                      :: T,Sigma,Omega
      real(kind = x_precision)                                 :: coeff_a=0.,coeff_b=0.,coeff_c=0.
      real(kind = x_precision),intent(out)                     :: H     
      real(kind = x_precision),intent(out)                     :: rho 
      real(kind = x_precision),intent(out)                     :: cs    
      real(kind = x_precision),intent(out)                     :: nu  
  
negative.code:     dichotomy             = (Smin+Smax)/2.
      j = 0
      call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs)
      call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_max,Sigma_0, Omega_0,rs)
      call variables(T, dichotomy, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
           K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_dichotomy,Sigma_0, Omega_0,rs)
      write(*,*)'fmin = ',f_min
      write(*,*)'fmax = ',f_max
      !call display_variables(T,Omega,r, sigma, H, rho, cs, nu, Q_plus, Q_minus,&
       !    K_ff, K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f)
      if ( f_max * f_min .gt. 0.) then
         write(*,*)'This function image does not switch its sign in this particular interval.'
      endif
      if( f_max * f_min .lt. 0.) then
         iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. 1000)
            if(f_min * f_dichotomy .lt. 0.) then
               Smax         = dichotomy
               call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
                    tau_eff, P_rad, P_gaz,E_ff,Fz,f_max,Sigma_0, Omega_0,rs)
            else
               Smin         = dichotomy
               call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
                    tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs)
            endif
            dichotomy       = (Smin+Smax)/2.0d0
            call variables(T, dichotomy, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
                 K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_dichotomy,Sigma_0, Omega_0,rs)
            j               = j + 1
         end do iteration
      endif
  