anchor.diff:      implicit none
       integer::i = 0
       real(kind = x_precision)                               :: temp = 0.0d0
  -    real(kind = x_precision), parameter                    :: t_min = 1.0d-3
  -    real(kind = x_precision), parameter                    :: t_max = 1.0d0
  +    real(kind = x_precision), parameter                    :: t_min = 1.0d-6
  +    real(kind = x_precision), parameter                    :: t_max = 2.0d-6
       integer,                  parameter                    :: nb_it = 10
       real(kind = x_precision)                               :: sigma = 0.0d0
  -    real(kind = x_precision)                               :: Smin = 1.0d-3
  -    real(kind = x_precision)                               :: Smax = 1.0d0
  -    real(kind = x_precision)                               :: eps = 1.0d-5
  +    real(kind = x_precision)                               :: Smin = 0.0d-5
  +    real(kind = x_precision)                               :: Smax = 0.0d-2
  +    real(kind = x_precision)                               :: eps = 1.0d-4
       real(kind = x_precision)                               :: omega = 0.0d0
  -    real(kind = x_precision)                               :: r = 0.0d0
  +    real(kind = x_precision)                               :: r = 5.0d0
       type(parameters)                                       :: param
       !----------------------------------------------------------------------
       
  
positive.code:   ! T-> Fixed variable
    !-------------------------------------------------------------------------
    dichotomy             = (Smin+Smax)/2.
  
    if ( f(T,Smin,omega) * f(T,Smax,omega) .gt. 0.) then
       write(*,*)'This function image does not switch its sign in this particular interval.'
    endif
    if( f(T,Smin,omega) * f(T,Smax,omega) .lt. 0.) then
       !iteration:do while ( dabs( f(T,dichotomy,omega) ) .ge. eps .and. j .lt. N)
       iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. N)
          if( f(T,Smin,omega) * f(T,dichotomy,omega) .lt. 0.) then
             Smax         = dichotomy
             write(*,*)'smin=',smin,' smax=', smax
  
negative.code:   else
       optical_depth     = 0
    end if
  write(*,*)'tau',tau_eff
    select case(optical_depth)
    case(1)
       Fz = 4._x_precision * c**2 * T**4/(27. * sqrt(3._x_precision) * (K_ff + K_e) * Sigma * Sigma_0)
  