c-------------------------------------------------------------------- 
c                        Fit halo profiles,
c                        use simple minimization
c-------------------------------------------------------------------- 
      PARAMETER (Nmax = 200000) ! 
      PARAMETER (Nbinmax = 85) !
      PARAMETER (Nqual   = 13) ! 
      COMMON  /HALOS/hals(Nqual,Nmax),Nhbin(Nmax),
     &       Prof(7,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)
      Character  Header*100

      Open(1,file='CatalogA.999.DAT')
      Open(2,file='CatFitAf.999.dat')
      Open(3,file='CatOut.999.dat')

      Lines =14    ! number of lines in the header preceeding radial bins
      lbins =85    ! number of bins along radius
      Lskip =4     ! number of header lines between bins and catalog

      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
         write (2,'(a)')  Header
         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)
         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.1.e12)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.
         Else
            CALL HaloFitExp(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc,
     &                      Rs(ih),VVmax(ih),r1(ih),rho1(ih),x2c(ih),
     &                      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.e0)
      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))
     &    ,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,g11.4,f7.2,f7.4)

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

 50   format(f8.3,i7,2g11.4,2f9.2)
            
             
       stop
      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-------------------------------------------------------------------- 
      PARAMETER (Nbinmax = 85) !
      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-------------------------------------------------------------------- 
      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-------------------------------------------------------------------- 
      PARAMETER (Nbinmax = 85) !
      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.
      rho1  = 1.022e3*(Umax/r1)**2 ! initial guess for rho_1
      dRho  = rho1/10.
      x2c   = 5.   
      a2c   = 1./10. 
      dx2   = x2c/10.
      da2   = a2c/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
         a2c_old  = a2c
         x2c_old  = x2c
         e_old    = error
c         write (*,*)  ' ==== iteration=== ',iter,error
      CALL LocalFit(Radden,Radout,Dens,istart,imax,Nbin,
     &        r1,rho1,a2c,x2c,dR,dRho,dx2,da2,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+
     &             (a2c/a2c_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
           a2c = a2c_old 
           x2c = x2c_old 
         dR =dR/1.21
         dRho =dRho/1.21
         dx2 =dx2/1.21
         da2 =da2/1.15
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
         da2 =min(da2/1.25,a2c/20.)
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,a2c,0.05)
      Rout =  Radden(imax)
      rhs  =  RhoExp(Rout,r1,rho1,x2c,a2c)/4.
      CALL FindRs(r_s2,r1,rho1,x2c,a2c,0.5,rhs,Rout)
      r_s =(r_s+r_s2)/2.
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,a2c,x2c,dR,dRho,dx2,da2,error,iFlag)
      PARAMETER (Nbinmax = 85) !
      PARAMETER (nStepR  = 6)   ! number of bins in radius 
      PARAMETER (nStepRho= 6)   ! number of bins in Rho
      PARAMETER (nStepa2 = 6)   ! number of bins in a2
      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.e-2)
      a2min = max(a2c - da2*nStepa2/2.,a2c/10.)

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," amin/max  =",4g11.4)')
c     &          a2min,a2min +nStepa2*da2,a2c 
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

      Do ia2  = 0,nStepa2
         a2   = a2min +ia2*da2
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,a2)/
     &                         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 
            a_2   = a2
            x_2   = x2
            ir_1  = ir
            irho_1= irho
            ix_1  = ix2
            ia_1  = ia2
        EndIf 
      EndDo 
      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.
      if(ia_1.eq.0.or.ia_1.eq.nStepa2)on_boundary = .true.
      iFlag = 0
      If(on_boundary)iFlag =1
      Rc  = r_1
      Rhoc= rho_1
      a2c = a_2
      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                         Exp-fit density profile
      FUNCTION RhoExp(r,r1,rho_1,x2,a2)
c                 input:  
c                          Vvmax = estimate of Vmax  
c-------------------------------------------------------------------- 
      x = r/r1
      RhoExp = rho_1/x*(exp(-x)+a2*exp(-x/x2))+1.

      Return
      End
c-------------------------------------------------------------------- 
c                         Exp-fit mass profile: dimensionless
c                         F =1-(1+x)exp(-x)
      FUNCTION FExp(x)
c-------------------------------------------------------------------- 
      real*8 y,Yexp
      y= x
      YExp = 1.d0-(1.d0+y)*exp(-y)
      FExp = YExp
      Return
      End
c-------------------------------------------------------------------- 
c                         Exp-fit velocity profile
      FUNCTION VcExp(r,r1,rho_1,x2,a2)
c                 output:  Vcirc in km/s
c                          
c-------------------------------------------------------------------- 
      x = r/r1
      VcExp = 6.698e-2*r1*sqrt(rho_1/x*(FExp(x)+a2*x2**2*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,a2,dr)
c                 output:  Vcirc in km/s
c                          
c-------------------------------------------------------------------- 

      Vmax =0.
      iter =0
      r    = r_1
      Vv0  = VcExp(r,r_1,rho_1,x2,a2)
      Vv1  = VcExp(r+dr,r_1,rho_1,x2,a2)
      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,a2)
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
c         write (*,*)  ' Vmax=',Vmax,' r=',r
c         r0 = r-dr/2.
c         V0= VcExp(r0,r_1,rho_1,x2,a2)  
c         r1 = r-dr
c         V1= VcExp(r1,r_1,rho_1,x2,a2)        
c         write (*,*)  ' V0=',V0,' r0=',r0      
c         write (*,*)  ' V1=',V1,' r1=',r1
      
      Vmax =0.945*Vcmax
      V2   = Vmax
      iter =0
      r    =0.
c      write (*,*)  ' V(rs)=',V2
 20   iter =iter +1    ! Find r_s: Vc(rs)=0.945Vmax
      r = r +dr
      Vv= VcExp(r,r_1,rho_1,x2,a2)
         If(Vv.lt.Vmax.and.r.lt.1.e4)Then
            GoTo 20
         EndIf 
         r_s =r
         r   =rmax
      Return
      End
c-------------------------------------------------------------------- 
c               Using density profile find Vmax and effective r_s for Exp-fit velocity profile
      SUBROUTINE FindRs(r_s,r_1,rho_1,x2,a2,dr,rhs,Rout)
c                 output:  Vcirc in km/s, r_s (kpc)
c               solve eq: rho(rs) = rho(Rout)/4*Rout/rs*(1+Rout/rs)**2               
c-------------------------------------------------------------------- 

      Vmax =0.
      iter =0
c      write (*,*)  ' Rout=',Rout,' rhs=',rhs
      r    = r_1/2.
 10   rho  = RhoExp(r,r_1,rho_1,x2,a2)
      rho_a= rhs*Rout/r*(1.+Rout/r)**2
c      write (*,*)  '     rho=',rho,' rho_a=',rho_a,r
      if(rho_a.lt.rho)Then
         r =r/2.
         If(r.lt.1.e-2)Then
            write (*,*)   'FindRs: initial guess failed'
            r_s = r_1/2
            return
         EndIf 
         goto 10
      EndIf 

 20   iter =iter +1    ! Find Vcmax using the approximation
      r = r +dr
      rho  = RhoExp(r,r_1,rho_1,x2,a2)
      rho_a= rhs*Rout/r*(1.+Rout/r)**2
c      write (*,*)  '     rho=',rho,' rho_a=',rho_a,r
      if(rho_a.gt.rho.and.r.lt.1.e3)goto20


         r_s =r
      Return
      End
