anchor.diff:      cs                   = Omega * H
       nu                   = params%alpha * cs * H
       K_ff                 = 6.13d22 * rho_0 * rho * (T_0 * T)**(-3.5_x_precision)
  -    K_e                  =  params%kappa_e
  +    K_e                  = params%kappa_e
       E_ff                 = 6.22d20 * (rho_0 * rho)**2 * sqrt(T_0 * T)
       tau_eff              = 0.5_x_precision * sqrt(K_e * K_ff) * Sigma * Sigma_0
   
  
positive.code: 
           write(number_of_cell,'(I5.5)') k
           fid = 20 + k
  
  !********For the optically thin medium********
          do i          = 1, nb_it
            Smin        = 1d-1
            Smax        = 1d5
  
            temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
  
            optical_depth = 0
  
            sigma       = dichotomy(Smin, Smax, eps, temp, omega, sigma_0, Omega_0,rs, T_0, rho_0, optical_depth)
  
  
            call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
                K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f,Sigma_0, Omega_0,rs, T_0, rho_0, optical_depth)
  
          !  call display_variables(temp,Omega,r, sigma, H, rho, cs, nu, Q_plus, Q_minus,&
          !      K_ff, K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f)
  
            temp_real_0(i)   = log10(temp * T_0)
            sigma_real_0(i)  = log10(sigma * sigma_0)
  
            optical_depth = 1
  
            call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
                K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f,Sigma_0, Omega_0,rs, T_0, rho_0, optical_depth)
  
            temp_real_1(i)   = log10(temp * T_0)
            sigma_real_1(i)  = log10(sigma * sigma_0)
  
         enddo
  
  
      !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( sigma_real_1(i) - sigma_real_0(j) ) .gt. eps2)
               write(fid,'(1p,E12.6,4x,1p,E12.6,4x,1p,E12.6)')sigma_real_1(i),temp_real_1(i)
               l = j
            enddo
         enddo
      enddo
  
      do i          = l, 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
  
    end subroutine curve
  
  
    !-------------------------------------------------------------------------
    !Subroutine for resolving a quadratic equation as long as the solutions are
    !a set of real numbers
  
negative.code:           j               = j + 1
  
         end do iteration
           dichotomy = S_center
  
      endif
  
  
    end function dichotomy
  
  