!------------------------------------------------------------
!                Different statistics for Einasto
!                dm density profile
!                      
!------------------------------------------------------------
Module Pars
Real*8, parameter :: rhos = 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
                                     G      = 6.672d-8,            &    ! Grav constant
                                     alpha = 0.170d0             ! slope for einasto profile
End Module Pars
!----------------------------------------------------------
Program Einasto
  use Pars
  Real*8 r,aMass,aM
  Do i=0,10
     r =10.**(i/2.+1)
     aM = aMass(r)
     write (*,*) r,aM
  EndDo

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 aMass(r)
!              mass profile
  use Pars
  Real*8 x,r,Faux,ANS,er,A
  External Faux

     A = 1.d-7
     x = max(r/rs,A)
     write(*,*) r,x,A
     Call DGAUS8 (Faux, A, x, ER, ANS, IERR)
     aMass = 4.*3.14159d0*1.d9*rhos*rs**3*ANS

   End Function aMass