anchor.diff:    
   contains
     !-------------------------------------------------------------------------
  -  ! subroutine to find the result of equation Q+ + Q- = 0
  -  ! output : datafile with temperature and surface density
  -  ! Range of surface density [Smin, Smax] and temperature [t_min, t_max]
  -  ! values can be changed
  +  ! Subroutine in order to find the result of the equation Q+ - Q- = 0
  +  ! output : datafile with the temperature T and the surface density Sigma
  +  ! Range of surface density [Smin, Smax] and temperature [T_min, T_max]
  +  ! values should be changed in order to find more coherent results
     !-------------------------------------------------------------------------
     subroutine curve()  
       implicit none
  -    integer::i
  -    real(kind = x_precision)                               :: temp
  -    real(kind = x_precision), parameter                    :: t_min = 1.0d7
  -    real(kind = x_precision), parameter                    :: t_max = 1.0d10
  -    integer,                  parameter                    :: nb_it = 3
  +    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
  +    integer,                  parameter                    :: nb_it = 10
       real(kind = x_precision)                               :: sigma = 0.0d0
  -    real(kind = x_precision)                               :: Smin = 1.0d7
  -    real(kind = x_precision)                               :: Smax = 1.0d10
  -    real(kind = x_precision)                               :: eps = 1.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)                               :: omega = 0.0d0
       real(kind = x_precision)                               :: r = 0.0d0
       type(parameters)                                       :: param
  
positive.code:     do i           = 1 , nb_it
         temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
         sigma       = dichotomy(Smin, Smax, eps, temp, omega)
         
         write(18,*)temp," ",sigma
      enddo
  
negative.code:      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
          else
             Smin         = dichotomy
             write(*,*)'smin=',smin,' smax=', smax
          endif
          dichotomy       = (Smin+Smax)/2.0d0
          j               = j + 1
  