anchor.diff:  
           ! S found with the dichotomy approach
           sigma_t_thick(i) = dichotomy(k, i, optical_depth)
  +        ! sigma_t_thick(i) = secant(k, i, optical_depth)
   
           ! Optical thin case (tau < 1)
           optical_depth = 0
   
           ! S found with the dichotomy approach
           sigma_t_thin(i)  = dichotomy(k, i, optical_depth)
  +        ! sigma_t_thin(i) = secant(k, i, optical_depth)
  +
         enddo
   
         ! For the first critical point (the one at the right of the S shape)
  
positive.code: 
      if ( f_max * f_min > 0.) then
  
        secante = 0
  
      else
  
        do while (dabs(Smax - Smin) >= eps .and. j < max_it)
  
           if (j == 0) then
              Smin = S_min
  
           else
              call variables(k, T, Smin, f_min, optical_depth, tau_eff)
              write(*,*)Smin
  
              Smin = Smin - (Smax - Smin)/(f_max - f_min) * f_min
           end if
  
          j = j + 1
  
        end do
  
       secante  = (Smin + Smax)/2._x_precision
  
      endif
  
    end function secante
  
  
  end module mod_s_curve
  
negative.code: 
  
    !-------------------------------------------------------------------------
    ! Secant method in order to determine the change of sign in a given
    ! interval [Smin,Smax] with an epsilon precision
    !-------------------------------------------------------------------------
    real(x_precision) function secant(k, i, optical_depth)
      implicit none
  
      integer, intent(in) :: k ! Position in the disk (r/Omega)
      integer, intent(in) :: i ! Iteration counter (Temperature)
      integer, intent(in) :: optical_depth ! Indicator for the optical thickness
      !-------------------------------------------------------------------------
      integer             :: j
      real(x_precision)   :: T ! Temperature used for the secant method
      real(x_precision)   :: tau_eff
      real(x_precision)   :: Smin     ! Lowest point
      real(x_precision)   :: Smax     ! Highest point
      real(x_precision)   :: f_min    ! Q+ − Q− at the lowest point
      real(x_precision)   :: f_max    ! Q+ − Q− at the highest point
      !-------------------------------------------------------------------------
  
      T = temperature(i)
  