c-------------------------------------------------------------------- 
c                        Fit halo profiles,
c                        use simple minimization
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (Nmax = 500000) ! 
      PARAMETER (Nbinmax = 120) !
      PARAMETER (Nqual   = 13) ! 
      COMMON  /HALOS/hals(Nqual,Nmax),Nhbin(Nmax),
     &       Prof(7,Nbinmax,Nmax),Nsmall(Nbinmax,Nmax)
      COMMON  /Binn/ Radout(Nbinmax)
      COMMON /APPROX/rsh(Nmax),Vmaxh(Nmax)
      DIMENSION Nin(Nbinmax),Radden(Nbinmax),Vc(Nbinmax),Dens(Nbinmax)
      DIMENSION VVmax(Nmax),Rs(Nmax),r1(Nmax),rho1(Nmax),x2c(Nmax),
     &               a2c(Nmax),error(Nmax),slope(Nmax)
      Character  Header*100

      Open(1,file='Catalog.unbound.DAT')
      Open(2,file='CatFit.unbound.dat')
      Open(3,file='CatOut.unbound.dat')
c      Open(1,file='CatalogB.486.DAT')
c      Open(2,file='CatFitB.486.dat')
c      Open(3,file='CatOut.486.dat')
c      Open(1,file='CatalogA.9999.DAT')
c      Open(2,file='CatFitA.9999.dat')
c      Open(3,file='CatOut.9999.dat')

      Lines =14    ! number of lines in the header preceeding radial bins
      lbins =100    ! number of bins along radius
      Lskip =5     ! number of header lines between bins and catalog
      aMfit = 1.e12   ! limit for mass for fitting NFW or Sersic profiles
      Do i=1,Lines
         read(1,'(a)') Header
         write (*,'(a)')  Header
         write (2,'(a)')  Header
      EndDo
      read(1,*) (Radout(i),i=1,lbins)
      write(2,'(10f8.2)') (Radout(i),i=1,lbins)
      write (*,'(10f8.2)') (Radout(i),i=1,lbins)
      Do i=1,Lskip
         read(1,'(a)') Header
         if(i.lt.4)write (2,'(a)')  Header
         if(i.eq.4)write (2,'(a,a)')  Header,
     &          '    C   Nbin   error Slope  nLarge'
         write (*,'(a)')  Header
      EndDo

      Nh =0
 100  Nh =Nh+1
      read(1,*,err=400,end=400) (hals(i,Nh),i=1,Nqual),Nhbin(Nh)
c         write(*,'(" Halo=",i8," coords=",3f8.3," bins=",i3)')
c     &              Nh,(hals(i,Nh),i=1,3),Nhbin(Nh)
         do i=1,Nhbin(Nh)
            read(1,*) (Prof(j,i,Nh),j=1,7),Nsmall(i,Nh)
         EndDo 
      If(Nh.lt.Nmax)goto100
c      If(Nh.lt.1000)goto100
      Nh =Nh +1
 400  Nh =Nh -1

      write (*,*)  ' Number of halos read =',Nh
      factor = 8.23e10
      Nchunk = max(Nh/20,100)
C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (ih,Umax,Rumax,Rvir,Nbin,Radden)
C$OMP+PRIVATE (Nin,Vc,Dens)
C$OMP+SCHEDULE(DYNAMIC,Nchunk)
      Do ih=1,Nh
         If(mod(ih,Nchunk/2).eq.0)write(*,*) ' halo=',ih
         Umax = hals(10,ih)
         Rumax= hals(12,ih)
         Rvir = hals(8,ih)
         Nbin = Nhbin(ih)
         do i=1,Nbin
            Radden(i) = Prof(1,i,ih)
            Nin(i)    = Prof(2,i,ih)
            Vc(i)     = Prof(6,i,ih)
            Dens(i)   = Prof(7,i,ih)/factor
         EndDo 
         if(hals(7,ih).lt.aMfit)Then
            CALL HaloFit(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc,
     &                  Rs(ih),VVmax(ih),Nbin,error(ih))
            r1(ih) = 0.
            rho1(ih) =0.
            x2c(ih) =0.
            a2c(ih) =0.
            slope(ih) =1.
         Else
            CALL HaloFitSER(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc,
     &                  Rs(ih),slope(ih),VVmax(ih),Nbin,error(ih))
c            CALL HaloFitExp(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc,
c     &                      Rs(ih),VVmax(ih),r1(ih),rho1(ih),x2c(ih),
c     &                      a2c(ih),Nbin,error(ih))
         EndIf
      EndDo 

      Do ih=1,Nh
         Rvir = hals(8,ih)
         Conc = Rvir/Rs(ih)
         error(ih) =min(error(ih),1.d0)
         Do i=1,Nhbin(ih)
            If(Prof(1,i,ih).ge.Rvir)Then
               nn =INT(Prof(2,i,ih))
               nL =INT(Prof(2,i,ih))-Nsmall(i,ih)
               goto 18
            EndIf 
         EndDo 
 18      write(2,20) (hals(i,ih),i=1,9),
     &    VVmax(ih),INT(hals(11,ih)),Rs(ih)*2.15,Conc,Nhbin(ih),
     &    sqrt(error(ih)),1./slope(ih),nL
c     &    ,r1(ih),rho1(ih),x2c(ih),a2c(ih)
 20   format(3f8.3,3F8.1,g11.3,f8.2,2f8.1,i9,f8.1,f8.2,i4,f8.3
     &       f7.3,i8)

c          If((Conc.lt.4..or.Conc.gt.25.).and.hals(7,ih).gt.1.e12)Then
c             write(3,'(" Halo=",i8," coords=",3f8.3," M=",g11.4,
c     &                 " C=",f8.3," R1,Rs,Rvir=",3f8.3)')
c     &                  ih,(hals(i,ih),i=1,3),hals(7,ih),
c     &                  Conc,r1(ih),Rs(ih),hals(8,ih)
c             write (3,'(5x,"r      Nin   Dens     Dens_4exp",
c     &              6x,"Vc  Vc_4exp")')
c        do i=1,Nhbin(ih)
c             Radd = Prof(1,i,ih)
c            Ni    = Prof(2,i,ih)
c            Vcc   = Prof(6,i,ih)
c            Den   = Prof(7,i,ih)/factor
c          
c         write(3,50) Radd,Ni,Den,
c     &             RhoExp(Radd,r1(ih),rho1(ih),x2c(ih),a2c(ih)),
c     &      Vcc,VcExp(Radd,r1(ih),rho1(ih),x2c(ih),a2c(ih))
c        EndDo 
c        EndIf
      enddo

 50   format(f8.3,i7,2g11.4,2f9.2)
            
             
       stop
      end
c-------------------------------------------------------------------- 
c                       Fit halo profile
      SUBROUTINE HaloFitSer(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc,
     &               Rs,slope,VVmax,Nbin,error)
c                 output:  Rs    = estimate of NFW core radius,
c                          Vvmax = estimate of Vmax
c                          slope = 1/n slope for Sersic fit  
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (Nbinmax = 120) !
      PARAMETER (facV    = 2.0) ! max variation for Vmax
      PARAMETER (facR    = 3.0)  ! max variation for rs
c      PARAMETER (facRvir = 2.0)  ! max radius for fit in units of Rvir
      PARAMETER (facRvir = 0.5)  ! max radius for fit in units of Rvir
      PARAMETER (nStepS  = 50)   ! number of bins in slope 
      PARAMETER (nStepR  = 50)   ! number of bins in radius 
      PARAMETER (nStepV  = 50)   ! number of bins in Vmax 
      DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin),
     &          Radout(Nbinmax)
      DATA iFlag/0/
      SAVE iFlag

      rs    = Rumax/2.15       ! initial guess for rs
      VVmax = Umax
      Vcmax = Umax
      Vmin  = Umax/facV
      dV    = Umax/nStepV*(facV -1./facV) ! step in Vmax
      Rmin  = rs/facR
      dR    = rs/nStepR*(facR -1./facR)   ! step in rs
      Smin  = 0.08 
      Smax  = 0.350
      dS    = (Smax-Smin)/nStepS          ! step in slope
      iFlag = iFlag +1

      error = 1.e10
      istart= 1
      Do istart =1,Nbin
         If(Nin(istart).gt.200)goto20
      EndDo 
         write (*,*)'not enough particles(first,last=)',Nin(1),Nin(Nbin)
      return           ! not enough particles; get back
 20   Rmax = min(Rvir/facRvir,Radout(Nbin))  ! max radius for the fit
      do i=Nbin,1,-1
         if(Radout(i).lt.Rmax)goto30
      EndDo 
 30   imax = i

      if(imax.le.istart)Then
         write (*,*)'not enough particles. Bins=',istart,imax
         return
      EndIf 
      Do is =0,nStepS
         slope  = Smin +is*dS
      Do iv =0,nStepV
         Vmaxx = Vmin +iv*dV
      Do ir =0,nStepR
         rss   = Rmin +ir*dR
         rho_0 = 1.022e3*(Vmaxx/rss)**2
         ee =0.
         ww =0.
         do i=istart,imax
          w  = (Radout(i)-Radout(i-1))/Radden(i) 
          ee=ee+w*(log10(RhoSER(Radden(i),rho_0,rss,slope)/Dens(i)))**2
          ww = ww + w
         EndDo
c         ee =ee/(imax-istart)
         ee =ee/ww
         If(ee.lt.error)then
            error = ee
            Rs    = rss
            Slp   = slope
            VVmax = Vmaxx
            ivv   = iv
            irr   = ir
            rho   = rho_0
c          write (*,'(10x,"V=",f8.3," r/slope=",2f8.3," err=",3g11.3)')
c     &          Vmaxx,rss,Slp,ee,rho_0
        EndIf 
      EndDo 
      EndDo 
      EndDo 
      write (*,'(10x,"V=",f8.3," r/Slp=",2f8.3," err=",7g11.3)')
     &          VVmax,Rs,1./Slp,error,rho,Rvir,Rvir/Rs
c      write(*,*) ' after ----'
c      do i=1,Nbin
c         x  = Radden(i)/rs
c         F  = log(1.+x)-x/(1.+x)
c         write(*,50) i,Radden(i),Dens(i),rho/x/(1.+x)**2,
c     &      Vc(i),VVmax*sqrt(4.62*F/x)
c      enddo
      slope = Slp
 50   format(i3,' r=',f8.3,' dens=',2g11.4,' v=',3f8.2)
      rho_0 = 1.022e3*(VVmax/Rs)**2
      Do i=1,Nbin
        write (3,'(4g11.4,i9)')Radout(i),Dens(i),Vc(i),
     &       RhoSER(Radden(i),rho_0,Rs,slope),Nin(i)
      EndDo 

c      if(iFlag.ge.10)stop
         Return
         End
c-------------------------------------------------------------------- 
c                         Sersic density profile
      FUNCTION RhoSER(r,rho_0,rs,slope)
c    
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
         x = r/rs
         RhoSER = rho_0/exp(2./slope*(x**slope-1.))+1.
      Return
      End
c-------------------------------------------------------------------- 
c                       Fit halo profile
      SUBROUTINE HaloFit(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc,
     &                  Rs,VVmax,Nbin,error)
c                 output:  Rs    = estimate of NFW core radius,
c                          Vvmax = estimate of Vmax  
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (Nbinmax = 120) !
      PARAMETER (facV    = 1.25) ! max variation for Vmax
      PARAMETER (facR    = 3.0)  ! max variation for rs
c      PARAMETER (facRvir = 2.0)  ! max radius for fit in units of Rvir
      PARAMETER (facRvir = 1.0)  ! max radius for fit in units of Rvir
      PARAMETER (nStepR  = 50)   ! number of bins in radius 
      PARAMETER (nStepV  = 20)   ! number of bins in Vmax 
      DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin),
     &          Radout(Nbinmax)

      rs    = Rumax/2.15       ! initial guess for rs
      VVmax = Umax
      Vcmax = Umax
      Vmin  = Umax/facV
      dV    = Umax/nStepV*(facV -1./facV) ! step in Vmax
      Rmin  = rs/facR
      dR    = rs/nStepR*(facR -1./facR)   ! step in rs

      error = 1.e10
      istart= 1
      Do istart =1,Nbin
         If(Nin(istart).gt.200)goto20
      EndDo 
      return           ! not enough particles; get back
 20   Rmax = min(Rvir/facRvir,Radout(Nbin))  ! max radius for the fit
      do i=Nbin,1,-1
         if(Radout(i).lt.Rmax)goto30
      EndDo 
 30   imax = i

      if(imax.le.istart)Then
c         write (*,*)'not enough particles. Bins=',istart,imax
         return
      EndIf 
      Do iv =0,nStepV
         Vmaxx = Vmin +iv*dV
      Do ir =0,nStepR
         rss   = Rmin +ir*dR
         rho_0 = 1.022e3*(Vmaxx/rss)**2
         ee =0.
         ww =0.
         do i=istart,imax
            w  = (Radout(i)-Radout(i-1))/Radden(i) 
            ee = ee + w*(log10(RhoNFW(Radden(i),rho_0,rss)/Dens(i)))**2
            ww = ww + w
         EndDo
c         ee =ee/(imax-istart)
         ee =ee/ww
         If(ee.lt.error)then
            error = ee
            Rs    = rss
            VVmax = Vmaxx
            ivv   = iv
            irr   = ir
            rho   = rho_0
c          write (*,'(10x,"V=",f8.3," r=",f8.3," err=",3g11.3)')
c     &          Vmaxx,rss,ee,rho_0
        EndIf 
      EndDo 
      EndDo 
c         write (*,'(10x,"V=",f8.3," r=",f8.3," err=",3g11.3)')
c     &          VVmax,Rs,error,rho
c      write(*,*) ' after ----'
c      do i=1,Nbin
c         x  = Radden(i)/rs
c         F  = log(1.+x)-x/(1.+x)
c         write(*,50) i,Radden(i),Dens(i),rho/x/(1.+x)**2,
c     &      Vc(i),VVmax*sqrt(4.62*F/x)
c      enddo
 50   format(i3,' r=',f8.3,' dens=',2g11.4,' v=',3f8.2)

         Return
         End
c-------------------------------------------------------------------- 
c                         NFW density profile
      FUNCTION RhoNFW(r,rho_0,rs)
c                 output:  Rs    = estimate of NFW core radius,
c                          Vvmax = estimate of Vmax  
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      x = r/rs
      RhoNFW = rho_0/x/(1.+x)**2

      Return
      End
c-------------------------------------------------------------------- 
c                       Double Exp-Fit for halo profile
      SUBROUTINE HaloFitExp(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc,
     &           Rs,VVmax,r1,rho1,x2c,a2c,Nbin,error)
c                 output:  Rs    = estimate of NFW core radius,
c                          Vvmax = estimate of Vmax  
c                 Uses 4 parameters:
c                         rho_1, r_1, a_2, x_2
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (Nbinmax = 120) !
      PARAMETER (facRho  = 10.) ! max variation for density
      PARAMETER (facR    = 10.0)  ! max variation for r1
c      PARAMETER (facRvir = 2.0)  ! max radius for fit in units of Rvir
      PARAMETER (facRvir = 0.5)  ! max radius for fit in units of Rvir
 
      DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin),
     &          Radout(Nbinmax)
      Logical   Converge

      rs    = Rumax/2.15       ! initial guess for rs
      VVmax = Umax
      r1    = min(max(Rumax/2.15,Radden(2)),Rvir/5.)       ! initial guess for r1
      dR    = r1/20.
c      rho1  = 1.022e3*(Umax/r1)**2 ! initial guess for rho_1
      rho1  = 200./(r1/Rvir)**2
      dRho  = rho1/10.
      x2c   = 1.5   
      dx2   = x2c/10.
      error = 100.
c      write (*,*)  
c      write(*,'(" Guess: R,rho,a2,x2=",4g11.3)')r1,rho1,a2c,x2c
c      write(*,'(" steps: R,rho,a2,x2=",4g11.3)')dR,dRho,da2,dx2

      istart= 1
      Do istart =1,Nbin
         If(Nin(istart).gt.100)goto20
      EndDo 
      return           ! not enough particles; get back
 20   Rmax = min(Rvir/facRvir,Radout(Nbin))  ! max radius for the fit
      do i=Nbin,1,-1
         if(Radout(i).lt.Rmax)goto30
      EndDo 
 30   imax = i
c      write (*,*)  ' max bin=',imax,' Rout=',Radout(imax),Rmax
      if(imax.le.istart.or.Nin(imax).lt.1000)Then
c         write (*,*)'not enough particles. Bins=',istart,imax
         return
      EndIf 

      Do iter =1,100
         r1_old   = r1
         rho1_old = rho1
         x2c_old  = x2c
         e_old    = error
c         write (*,*)  ' ==== iteration=== ',iter,error
      CALL LocalFit(Radden,Radout,Dens,istart,imax,Nbin,
     &        r1,rho1,x2c,dR,dRho,dx2,error,iFlag)
c          write (*,'(10x,"r1=",f8.4," rho1=",7g11.4," err=",g11.4,i2)')
c     &          r1,rho1,a2c,x2c,dR,dRho,dx2,da2,error,iFlag
         quality = (r1/r1_old-1.)**2+(rho1/rho1_old-1.)**2+
     &             (x2c/x2c_old-1)**2
         Converge = (quality.lt.1.e-4).and.
     &              (abs(error/e_old-1.).lt.1.e-3).and.
     &              (iter.gt.10).and.(iFlag.eq.0)
        dRho =max(min(dRho,rho1/10.),rho1/30.) 
         If(Converge)Goto40
      If(error.ge.e_old.and.iter.gt.10)Then
           r1  = r1_old  
           rho1= rho1_old
           x2c = x2c_old 
         dR =dR/1.21
         dRho =dRho/1.21
         dx2 =dx2/1.21
c          write (*,*)  '        restore center, reduce steps'
      EndIf
      If(iFlag.eq.0)Then
         dR =min(dR/1.41,r1/20.)
         dRho =min(dRho/1.41,rho1/20.)
         dx2 =dx2/1.41
c         write (*,*)  '        changed steps'
c      Else
c         dR =dr/1.01
c         dRho =drho/1.01
      EndIf 
      EndDo 

 40   Rs =r_1
c       write (*,'(10x,"r1=",f8.3," rho1=",7g11.4," err=",g11.3,2i4)')
c     &          r1,rho1,a2c,x2c,dR,dRho,dx2,da2,error,iFlag,iter

      CALL  FindVcMax(rmax,Vcmax,Rvir,r_s,r1,rho1,x2c,0.05)
      Rout =  Radden(imax)
c      r_s =min(r_s,r_s2)
c      r_s =r_s2
c      If(Rvir/r_s.lt.4.)Then
c      write (*,'(" rmax,Vmax=",2f8.3," rs=",3g11.3)')rmax,Vcmax,
c     &                         r_s,Rout,Dens(imax)
c      do i=istart,imax
c         write(*,50) i,Radden(i),Dens(i),
c     &             RhoExp(Radden(i),r1,rho1,x2c,a2c),
c     &            rhs*Rout/Radden(i)*(1.+Rout/Radden(i))**2,
c     &      Vc(i),VcExp(Radden(i),r1,rho1,x2c,a2c)
c      enddo
c         EndIf
c 50   format(i3,' r= ',f8.3,' dens=',3g11.4,' v=',3f8.2)
      Rs    = r_s
      VVmax = Vcmax

      Return
      End
c-------------------------------------------------------------------- 
c                       Fit halo profile
      SUBROUTINE LocalFit(Radden,Radout,Dens,istart,imax,Nbin,
     &        Rc,Rhoc,x2c,dR,dRho,dx2,error,iFlag)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (Nbinmax = 120) !
      PARAMETER (nStepR  = 6)   ! number of bins in radius 
      PARAMETER (nStepRho= 6)   ! number of bins in Rho
      PARAMETER (nStepx2 = 6)   ! number of bins in x2
      DIMENSION Radden(Nbin),Dens(Nbin),Radout(Nbinmax)
      Logical on_boundary

      Rmin = max(Rc -dR*nStepR/2.,Rc/10.)
      Rhomin = max(Rhoc -dRho*nStepRho/2.,Rhoc/10.)
      x2min = max(x2c - dx2*nStepx2/2.,1.d-2)

c      write (*,'(5x," Rmin/max  =",4g11.4)')
c     &          Rmin,Rmin +nStepR*dR,Rc 
c      write (*,'(5x," Rhomin/max=",4g11.4)')
c     &          Rhomin,Rhomin +nStepRho*dRho,Rhoc 
c      write (*,'(5x," xmin/max  =",4g11.4)')
c     &          x2min,x2min +nStepx2*dx2,x2c 
      error =1.e+10
      Do ir =0,nStepR
         r11   = Rmin +ir*dR
         e_old = error
      Do irho =0,nStepRho
         rho11   = Rhomin +irho*dRho
      
      Do ix2  = 0,nStepx2
         x2   = x2min +ix2*dx2

c         write (*,'(4g11.3,4i4)') r11,rho11,x2,a2,ir,irho,ix2,ia2
         ee =0.
         ww =0.
         do i=istart,imax
            w  = (Radout(i)-Radout(i-1))/Radden(i) 
            ee = ee + w*(log10(RhoExp( Radden(i),r11,rho11,x2)/
     &                         Dens(i)))**2
            ww = ww + w
         EndDo
c         ee =ee/(imax-istart)
         ee =ee/ww
         If(ee.lt.error)then
            error = ee
            r_1   = r11
            rho_1 = rho11 
            x_2   = x2
            ir_1  = ir
            irho_1= irho
            ix_1  = ix2
        EndIf 
      EndDo 

      EndDo
c          write (*,'(10x,"r1=",f8.3," rho1=",3g11.3," err=",g11.3,4i4)')
c     &          r_1,rho_1,a_2,x_2,error,ir_1,irho_1,ix_1,ia_1
c         If(abs(error-e_old).lt.1.e-5)Goto100 
      EndDo 

 100  on_boundary = .false.
      if(ir_1.eq.0.or.ir_1.eq.nStepR)on_boundary = .true.
      if(irho_1.eq.0.or.irho_1.eq.nStepRho)on_boundary = .true.
      if(ix_1.eq.0.or.ix_1.eq.nStepx2)on_boundary = .true.
      iFlag = 0
      If(on_boundary)iFlag =1
      Rc  = r_1
      Rhoc= rho_1
      x2c = x_2
c          write (*,'(10x,"r1=",f8.3," rho1=",3g11.3," err=",g11.3,4i4)')
c     &          r_1,rho_1,a_2,x_2,error,ir_1,irho_1,ix_1,ia_1

         Return
         End

c-------------------------------------------------------------------- 
c                         LogNorm-fit for density profile
      FUNCTION RhoExp(r,r1,rho_1,x2)
c
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      x = r/r1
      RhoExp = rho_1/x**2*exp(-0.5*(x/x2)**2)+1.

      Return
      End
c-------------------------------------------------------------------- 
c                        dimensionless integrant for mass profile
      FUNCTION AxExp(x,x2)
c                 
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
         AxExp = exp(-0.5*(LOG(x)/x2)**2)
      Return
      End
c-------------------------------------------------------------------- 
c                         Exp-fit mass profile: dimensionless
c                        
      FUNCTION FExp(x,x2)
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8    INTG
      EXTERNAL  AxEXP
         xmin = min(1.d-4,x/10.d0)
         FExp = INTG(AxEXP,xmin,x,x2)
      Return
      End
c-------------------------------------------------------------------- 
c                         Exp-fit velocity profile
      FUNCTION VcExp(r,r1,rho_1,x2)
c                 output:  Vcirc in km/s
c                          
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
         x = r/r1
         VcExp = 6.698e-2*r1*sqrt(rho_1/x*FExp(x,x2))
      Return
      End
c-------------------------------------------------------------------- 
c                        Vmax and effective r_s for Exp-fit velocity profile
      SUBROUTINE FindVcMax(r,Vcmax,Rvir,r_s,r_1,rho_1,x2,dr)
c                 output:  Vcirc in km/s
c                          
c-------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      Vmax =0.
      iter =0
      r    = r_1
      Vv0  = VcExp(r,r_1,rho_1,x2)
      Vv1  = VcExp(r+dr,r_1,rho_1,x2)
      If(Vv1.gt.Vv0)Then
         r = r-dr
      Else
         r = 0
      EndIf
 10   iter =iter +1    ! Find Vcmax using the approximation
      r = r +dr
      Vv= VcExp(r,r_1,rho_1,x2)
c      write (*,*)  '         VV,r=',Vv,r,Vmax
         If(Vv.gt.Vmax.and.r.lt.1.e4.and.r.lt.Rvir)Then
            Vmax = Vv
            GoTo 10
         EndIf 

         Vcmax = Vmax
         rmax  = r
      Return
      End
C-------------------------------------------- Simpson integration
      REAL*8 FUNCTION INTG(FUNC,A,B,sig)
C---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (EPS=5.0d-4, JMAX=22)
      EXTERNAL FUNC
      OST=-1.D30
      OS= -1.D30
      ST =0.
      DO 11 J=1,JMAX
        CALL TRAPZD(FUNC,A,B,ST,J,sig)
        INTG=(4.0d0*ST-OST)/3.0d0
        IF (ABS(INTG-OS).Le.EPS*ABS(OS)) RETURN
        OS=INTG
        OST=ST
11    CONTINUE
      WRITE (*,*)'Integration did not converge'
      RETURN
      END
C----------------------------------------------
      SUBROUTINE TRAPZD(FUNCC,A,B,S,N,sig)
C---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
        SAVE IT
        EXTERNAL FUNCC
      IF (N.EQ.1) THEN
        S=0.5d0*(B-A)*(FUNCC(A,sig)+FUNCC(B,sig))
        IT=1
      ELSE
        TNM=IT
        DEL=(B-A)/TNM
        X=A+0.5D0*DEL
        SUM=0.0D0
        DO 11 J=1,IT
          SUM=SUM+FUNCC(X,sig)
          X=X+DEL
11      CONTINUE
        S=0.5D0*(S+(B-A)*SUM/TNM)
        IT=2*IT
      ENDIF
      RETURN
      END

