anchor.diff:      real(kind = x_precision),intent(out)                     :: P_rad 
       real(kind = x_precision),intent(out)                     :: P_gaz 
       real(kind = x_precision),intent(in)                      :: Sigma_0
  +    real(kind = x_precision),intent(in)                      :: Omega_0
  +    real(kind = x_precision),intent(in)                      :: rs
       real(kind = x_precision),intent(out)                     :: f 
       integer                                                  :: optical_depth = 0
       !------------------------------------------------------------------------
   
       call get_parameters(param)
      
  -    coeff_a              = (Omega**2 * state_0%Omega_0**2 * Sigma * Sigma_0)/2.
  +    coeff_a              = (Omega**2 * Omega_0**2 * Sigma * Sigma_0)/2.
       coeff_b              = (-1._x_precision/3._x_precision) * cst_rad*T**4 * state_0%T_0**4
       coeff_c              = (-1._x_precision * param%RTM * T * Sigma * Sigma_0)/2._x_precision
  -  
  +
       call quadratic(coeff_a,coeff_b,coeff_c,H)
  -  
       rho                  = Sigma/H
  -    cs                   = Omega * H
  +    !cs                   = Omega * H
  +    cs                   = Omega_0 * rs
       nu                   = param%alpha * cs * H
       P_rad                = T**4
       P_gaz                = rho * T
  
positive.code:     !------------------------------------------------------------------------
      call get_parameters(param)
      call display_parameters()
      call initial_variables(rs, rmin, Mdot_0, Sigma_0)
      call display_initial_variables(rs, rmin, Mdot_0, Sigma_0)
      !-------------------------------------------------------------------------
      ! Test for 1 value of r
      !-------------------------------------------------------------------------
      do j              = 1 , 1
         !do j              = 1 , n_cell
         !  r = (param%rmax-rmin)/(n_cell-1)*(j-1) + rmin
         r              = 6*G*param%M/(c**2)     
         omega          = G*param%M/r**3
         write(number_of_cell,'(I5.5)') j
         fid = 20 + j
         fname = 'Temperature_Sigma_'//trim(number_of_cell)//'.dat'
         
         open(fid,file   = fname, status='unknown',action='readwrite') 
         write(fid,*)"T              Sigma"
         do i           = 1 , nb_it
            Smin        = 1.0d-5
            Smax        = 1.0d-0
            temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
            sigma       = dichotomy(Smin, Smax, eps, temp, omega)
            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)
            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)
            
            write(fid,'(1p,E12.6,4x,1p,E12.6)')temp,sigma
         enddo
         close(fid)
  
negative.code:   ! Function in order to determine the change of sign in a given interval
    ! [Smin,Smax] with an epsilon precision
    !-------------------------------------------------------------------------
    real(kind=x_precision) function dichotomy(Smin, Smax, eps, T, omega, sigma_0, Omega_0,rs)
      use mod_read_parameters
      use mod_constants
      use mod_variables
  