anchor.diff:  
             sigma       = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
   
  -          call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
  -               K_e, tau_eff, E_ff,Fz,f, optical_depth)
  -
  -        !  call display_variables(temp,Omega,r, sigma, H, rho, cs, nu, Q_plus, Q_minus,&
  -        !      K_ff, K_e, tau_eff, E_ff,Fz,f)
  +          call variables(temp, sigma, omega, f, optical_depth)
   
             temp_t_0(i)   =  temp
             sigma_t_0(i)  =  sigma
  
positive.code:     coeff_c              = - params%RTM * T * Sigma * state_0%S_0 / (2._x_precision * state_0%H_0**2)
  
      call quadratic(coeff_a , coeff_b , coeff_c , H)
  
  
      rho                  = Sigma / H
      cs                   = Omega * H
  
negative.code:     !-------------------------------------------------------------------------
      S_center             = (Smin+Smax)/2.
      j = 0
      call variables(T, Smin, Omega, f_min, optical_depth)
  
      call variables(T, Smax, Omega, f_max, optical_depth)
  
    !   write(*,*)'fmin = ',f_min
    !   write(*,*)'fmax = ',f_max
  