anchor.diff:  
             call variables(temp, sigma, omega, f, optical_depth)
   
  -          temp_t_1(i)   =  temp
  -          sigma_t_1(i)  = sigma
  +          temp_t_thick(i)   =  temp
  +          sigma_t_thick(i)  = sigma
   
             !do j          = 1, nb_it
             !   write(fid_2,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_1(j),temp_real_1(j)
  
positive.code:             ! write(*,*)index_fcp
  
      ! do i = index_fcp, nb_it
      !   write(*,*)sigma_real_0(i), sigma_real_1(i)
      !   write(*,*)(dabs(sigma_real_1(i) - sigma_real_0(i+1)))
      ! end do
  
      do i = index_fcp, nb_it - 1
  
         if(dabs(sigma_real_1(i) - sigma_real_0(i+1)) .lt. eps2)then
            index_scp = i
            ! sigma_c_thin = (sigma_real_1(i) + sigma_real_0(i+1))/2._x_precision
            sigma_c_thin =  sigma_real_0(i+1)
  
            temp_c_thin = temp_real_1(i)
         end if
      enddo
  
  
negative.code:        enddo
  
  
         call first_critical_point(sigma_t_thick, temp_t_thick, index_fcp,sigma_c_thick, temp_c_thick, nb_it)
  
         call second_critical_point(sigma_t_thick, sigma_t_thin, temp_t_thick,&
           index_fcp, index_scp, sigma_c_thin, temp_c_thin, nb_it, eps2)
  
         call build_s_curve(sigma_t_thick, sigma_t_thin, temp_t_thick, nb_it, index_scp, sigma_real, temp_real)
  
         call display_critical_points(sigma_c_thin, temp_c_thin,sigma_c_thick, temp_c_thick, k)
  
  