anchor.diff:      real(kind = x_precision),intent(out)                     :: Q_minus
       real(kind = x_precision),intent(out)                     :: P_rad
       real(kind = x_precision),intent(out)                     :: P_gaz
  -
       real(kind = x_precision),intent(out)                     :: f
  -
       !------------------------------------------------------------------------
   
       coeff_a              = (Omega**2 * Omega_0**2 * Sigma * Sigma_0)/2._x_precision
  
positive.code:     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,T_0, rho_0, optical_depth)
  
      ! write(*,*)'fmin = ',f_min
      ! write(*,*)'fmax = ',f_max
  
      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. 10000)
  
  
  
negative.code:     real(kind = x_precision),intent(in)                      :: rho_0
      real(kind = x_precision),intent(in)                      :: rs
      integer, intent(in)                                      :: optical_depth
      real(kind = x_precision)                                 :: coeff_a=0.,coeff_b=0.,coeff_c=0.
  
      real(kind = x_precision),intent(out)                     :: H
  