anchor.diff:            temp_real_0(i)   = log10(temp * T_0)
             sigma_real_0(i)  = log10(sigma * sigma_0)
   
  -       enddo
  -
  -
  -       temp = 0.0d0
  -       sigma = 0.0d0
  -
  -!********For the optically thick medium********
  -        do i          = 1, nb_it
  -          Smin        = 1d-2
  -          Smax        = 1d51
  -          t_min       = 3.d-1
  -          t_max       = 2.49d0
  -
  -          temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
  -          sigma       = dichotomy(Smin, Smax, eps, temp, omega, sigma_0, Omega_0,rs, T_0, rho_0)
  -
             optical_depth = 1
   
             call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
  
positive.code: 
      !Saving all the data
      fname = 'Temperature_Sigma_'//trim(number_of_cell)//'_r.dat'
         open(fid,file  = fname, status='unknown',action='readwrite')
         !  write(fid,*)"T     &         Sigma"
         do i          = 1, nb_it
            do j          = 1, nb_it
               do while( abs( temp_real_1(i) - temp_real_0(j) ) .gt. eps .AND. &
                          abs( sigma_real_1(i) - sigma_real_0(j) ) .gt. eps) then
  
                  write(fid,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_1(i),temp_real_1(i)
  
                enddo
             enddo
          enddo
  
          do i          = j, nb_it
  
               write(fid,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_0(i),temp_real_0(i)
  
          enddo
         close(fid)
  
      enddo
  
  
negative.code: 
  
      call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs,T_0,rho_0, optical_depth)
  
      call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_max,Sigma_0, Omega_0,rs,T_0, rho_0, optical_depth)
  
      call variables(T, S_center, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
           K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_center,Sigma_0, Omega_0,rs,T_0, rho_0, optical_depth)
  
  
            if(f_min * f_center .gt. 0.) then
  