!------------------------------------------------------------
!                Different statistics for Einasto
!                dm density profile
!                      
!------------------------------------------------------------
Module Pars
    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
  Real*8 x,fMass,f,y,v2,Rho,alpha0,dalpha,app,app2,app3,C0,C1,C2, &
         err,err2,ee,eru2,Cc0,Cc1,Cc2
write(*,*) ' X   M(<x)    y       Vcirc             dens/<rho>'



eru  = 1.e10
eru2 = 1.e10
Do k1 = -2,2
   C0 = 3.0551-k1*0.0002d0
   write(*,*) '  --- C0=',C0,k1
Do k2 = -5,5
   C1 = 0.63050-k2*0.0002d0 
Do k3 = -5,5
   C2 = 0.3727-k3*0.0002d0 

err  = 0.
err2 = 0.   
alpha0 = 0.06 ; dalpha =0.01
alpha = alpha0 
Do while(alpha<0.3)
   
  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 = fMass(x)
        v2 = exp(2./alpha)*f/x
        If(f .ge. y)exit
     EndIf
  EndDo
     app2= C0*exp(-C1*alpha**C2)
     ee  = (app2/x-1.)**2
     err = err + ee
     err2 = max(err2,ee)
        alpha = alpha +dalpha
end Do
     If(err2.lt.eru2)Then
        eru2 = err2
        Cc0 = C0 ; Cc1 =C1 ; Cc2 =C2
      write(*,'(a,3f8.4,a,f8.4)')' Best values =',Cc0,Cc1,Cc2,' Max Error =',sqrt(eru2)*100.        
     EndIf
End Do
end Do
End Do
write(*,'(a,3f8.4)') ' Best values =',Cc0,Cc1,Cc2
write(*,'(a,3f8.4)') ' Max Error   =',sqrt(eru2)*100.

alpha0 = 0.06 ; dalpha =0.01
alpha = alpha0 
Do while(alpha<0.3)
   
  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 = fMass(x)
        v2 = exp(2./alpha)*f/x
        If(f .ge. y)exit
     EndIf
  EndDo
  app2= Cc0*exp(-Cc1*alpha**Cc2)
        write(*,'(1p,12g13.5)') x,f,y,v2,Rho(x),alpha,  &
                   100.*(app2/x-1.),app2
        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
   End Function fMass