!------------------------------------------------------------
!                Different statistics for Einasto
!                dm density profile
!                      
!------------------------------------------------------------
Module Pars
  IMPLICIT REAL*8 (A-H,O-Z)
    Real*8, parameter :: rhos = 1.0*1.73d-3,     &    ! Rho_s in unit of Msun/pc**3
                         rs      = 21.5,          &    ! Rs in kpc
                         h0     = 0.73,          &    ! Hubble constant
                         rho_cr = 275.5*h0**2,  & ! critical density in Msun/kpc**3
                         Omeg_m =0.27,             & ! matter
                         G      = 6.672d-8,            &    ! Grav constant
                         pi      = 3.141592653d0
Real*8                   alpha                     ! slope for einasto profile

End Module Pars
!----------------------------------------------------------
Program Einasto
  use Pars
  IMPLICIT REAL*8 (A-H,O-Z)
  Real*8 fMass

alpha0 = 0.10 ; dalpha =0.02
alpha = alpha0 
  Cvir  = 10.

Do while(alpha<0.22)
   
  Do i=0,50000
     x =1.5d0*10.**(i/5000.d0)
     If(x<5.)Then
        y = x**3/exp(2.d0/alpha*x**alpha)
        f = 4.*fMass(x)
        v2 = exp(2./alpha)*f/x
        If(f .ge. y)exit
     EndIf
  EndDo
  Xmax  = x
  Fvir  = exp(2./alpha)*fMass(Cvir)
  Vmax  = sqrt(exp(2./alpha)*fMass(Xmax)/Xmax *(Cvir/Fvir))
   write(*,'(/1p,6(a,g12.4))') ' ------- alpha = ',alpha, &
                               ' Cvir =',Cvir,    &
                               ' Fvir=',Fvir,' Xmax =',Xmax, &
                               ' Vmax/Vrir =',Vmax
   write(*,'(2a)') '   x         M(<x)     Vcirc/Vvir      dens/<rho>', &
          ' NFW: M(<x)     Vcirc/Vvir      dens/<rho>' 
  Do i=0,20
     x =0.1d0*10.**(i/5.d0)
     If(x<50.)Then
        f = exp(2./alpha)*fMass(x)
        v2 = sqrt(f/x *(Cvir/Fvir))
        rho = 0.25/exp(2.d0/alpha*(x**alpha-1))
        fnfw = fMnfw(x)
        v = sqrt(fnfw/x *(Cvir/fMnfw(Cvir)))
        rhonfw = 1./x/(1.+x)**2


       write(*,'(1p,12g13.5)') x,f,v2,rho,fnfw,v,rhonfw
    end If
  EndDo
        alpha = alpha +dalpha
end Do

End Program Einasto
!----------------------------------------------------------
Real*8 Function Rho(x)
!
  use Pars
  Real*8 x
     Rho = exp(-2.d0/alpha*(x**alpha -1.d0))
   End Function Rho

!----------------------------------------------------------
Real*8 Function Faux(x)
!              integrant for mass profile
  use Pars
  Real*8 x

    Faux = exp(-2.d0/alpha*(x**alpha ))*x**2
  End Function Faux
!----------------------------------------------------------
Real*8 Function fMass(x)
!              mass profile
  use Pars
  Real*8 x,r,Faux,ANS,er,A
  External Faux

     Call DGAUS8 (Faux, 0.d0, x, ER, ANS, IERR)
     fMass = ANS/4.
   End Function fMass
!----------------------------------------------------------
Real*8 Function fMnfw(x)
!              mass profile for NFW
  use Pars
  Real*8 x

   fMnfw = log(1.+x) -x/(1.+x)
  
 End Function fMnfw