!---------------------------------------------------
!
!               Get statistics for a table of Lambda
!
!---------------------------------------------------
Program Stat
    Character :: Line*120 
    Integer, Parameter :: Nbin = 500, M =500000, Nmax =8000000
    Real :: DD(M),Vc(Nmax),Rmc(Nmax),aMc(Nmax), &
                    aMtot(Nmax),aLambda(Nmax),dR(Nmax),Vratio(Nmax)
    Integer :: indR(M),Indx(-Nbin:Nbin)
    Real      :: aLdistr(-Nbin:Nbin)
    pi = 3.1415926

Open(1,file='CATALOGS/CatshortDistinctd5.416.DAT')
Open(2,file='ListLambdaDistinctStat.dat')
Open(3,file='LambdaDistinct.dat')

Do i=1,11
   read(1,'(a)') Line
   write(2,'(a)')Line
EndDo

Ncenter =0
Do 
   Read(1,*,iostat=iStat)i,x,y,z,vx,vy,vz,aM,rr,vrms1,vv,aM2,vrms2,VirRat,ddr,aJ,aL
   If(iStat/= 0)Exit
      Ncenter = Ncenter +1
      If(Ncenter > Nmax)Stop ' Not enough space. Increase Nmax'
      Vc(Ncenter) = vv*sqrt(1.+(20/vv)**2)
      aMc(Ncenter) = aM
      Rmc(Ncenter) = rr
      aMtot(Ncenter) = aM2
      Vratio(Ncenter) = VirRat
      aLambda(Ncenter) = aL
EndDo
   write(*,*) ' Read halos=',Ncenter

   dlog = 0.05
   write(2,'(2a,i7)')'Nhalos       <M>     L1          L10         Lmedian       L90         L99  ',' N=',Ncenter
   aMin = 1.e9
   iFlag = 0
Do ibin =1,Nbin
   icount =1
10   aMax =aMin*10.**(dlog)
   iC =0
   ss =0.
   Do i=1,Ncenter
      If(aMc(i)>=aMin.and.aMc(i)<aMax.and.Rmc(i)>0.01)Then
         iC =iC +1
         If(iC>M)Stop ' Not enough space for buffer. Increase M'
         DD(iC) = aLambda(i)    ! store redults for ranking
         ss = ss + aMc(i)
      EndIf
   EndDo
   If(iC>5.and.iFlag==0)iFlag=1
   If(iC < 250.and.aMin<5.e14.and.iFlag==1.and.icount<10)Then
      write (*,*) iC,dlog,aMin
      dlog = dlog*1.1
      icount = icount +1
      goto 10
   EndIf
      If(iC>10)Then
      CALL R_mrgref (DD,indR,iC)  ! rank values
      M1  = iC/100+1
      M10 = iC/10+1
      M90 = 9*iC/10+1
      M50 = iC/2+1
      M99 = 99*iC/100+1
      V1  = DD(indR(M1))
      V10 = DD(indR(M10))
      V50 = DD(indR(M50))
      V90 = DD(indR(M90))
      V99 = DD(indR(M99))

     write(*,'(a,i8,a,1p,6g12.5)') ' Selected halos=',iC,' Masses=',aMin,aMax,V50
      write(2,'(i8,1p,6g12.5)') iC,ss/iC,V1,V10,V50,V90,V99
      EndIf
      aMin =aMax
EndDo
! ---------------------------------------------------------- Distribution of Lambda
   dlog  =  0.010
   aMin  = 1.3e8*500.
   aMax = 5.e13
   iC       = 0
   Indx   = 0
   aLdistr =0.
   Do i=1,Ncenter
      If(aMc(i)>=aMin.and.aMc(i)<aMax.and.Rmc(i)>0.01)Then
         iC =iC +1
          j  = INT(log10(aLambda(i))/dlog+1000)-1000    ! store redults for ranking
          j = Max(Min(Nbin,j),-Nbin)
          Indx(j) = Indx(j) + 1
          aLdistr(j) = aLdistr(j) + aLambda(i)
      EndIf
   EndDo
   write(3,*) ' ----------  Total number of halos selected =',iC
   write(3,*) '                MassMin/max =',aMin,aMax
   write(3,*) '                dlog_10            =',dlog
   write(3,*) ' Lmabda_min <Lambda> dN/dlog_10/N   dN'  
   iFlag =0
   Do j=-Nbin,Nbin
      if(Indx(j)>0)Then
         aL =aLdistr(j)/Indx(j)
         dst = Indx(j)/dlog/iC
         write(3,'(3f10.5,2i7)')10.**((j)*dlog),aL,dst,Indx(j)
      EndIf
   EndDo
End Program Stat

Subroutine R_mrgref (XVALT, IRNGT,N)
!   Ranks array XVALT into index array IRNGT, using merge-sort
! __________________________________________________________
!   This version is not optimized for performance, and is thus
!   not as difficult to read as some other ones.
!   Michel Olagnon - April 2000
! __________________________________________________________
! _________________________________________________________
      Real, Dimension (N), Intent (In) :: XVALT
      Integer, Dimension (N), Intent (Out) :: IRNGT
! __________________________________________________________
!
      Integer, Dimension (:), Allocatable :: JWRKT
      Integer :: LMTNA, LMTNC
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      NVAL = Min (SIZE(XVALT), SIZE(IRNGT))
      If (NVAL <= 0) Then
         Return
      End If
!
!  Fill-in the index array, creating ordered couples
!
      Do IIND = 2, NVAL, 2
         If (XVALT(IIND-1) < XVALT(IIND)) Then
            IRNGT (IIND-1) = IIND - 1
            IRNGT (IIND) = IIND
         Else
            IRNGT (IIND-1) = IIND
            IRNGT (IIND) = IIND - 1
         End If
      End Do
      If (Modulo (NVAL, 2) /= 0) Then
         IRNGT (NVAL) = NVAL
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      Allocate (JWRKT(1:NVAL))
      LMTNC = 2
      LMTNA = 2
!
!  Iteration. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (LMTNA >= NVAL) Exit
         IWRKF = 0
         LMTNC = 2 * LMTNC
         IWRK = 0
!
!   Loop on merges of A and B into C
!
         Do
            IINDA = IWRKF
            IWRKD = IWRKF + 1
            IWRKF = IINDA + LMTNC
            JINDA = IINDA + LMTNA
            If (IWRKF >= NVAL) Then
               If (JINDA >= NVAL) Exit
               IWRKF = NVAL
            End If
            IINDB = JINDA
!
!   Shortcut for the case when the max of A is smaller
!   than the min of B (no need to do anything)
!
            If (XVALT(IRNGT(JINDA)) <= XVALT(IRNGT(JINDA+1))) Then
               IWRK = IWRKF
               Cycle
            End If
!
!  One steps in the C subset, that we create in the final rank array
!
            Do
               If (IWRK >= IWRKF) Then
!
!  Make a copy of the rank array for next iteration
!
                  IRNGT (IWRKD:IWRKF) = JWRKT (IWRKD:IWRKF)
                  Exit
               End If
!
               IWRK = IWRK + 1
!
!  We still have unprocessed values in both A and B
!
               If (IINDA < JINDA) Then
                  If (IINDB < IWRKF) Then
                     If (XVALT(IRNGT(IINDA+1)) > XVALT(IRNGT(IINDB+1))) &
                    & Then
                        IINDB = IINDB + 1
                        JWRKT (IWRK) = IRNGT (IINDB)
                     Else
                        IINDA = IINDA + 1
                        JWRKT (IWRK) = IRNGT (IINDA)
                     End If
                  Else
!
!  Only A still with unprocessed values
!
                     IINDA = IINDA + 1
                     JWRKT (IWRK) = IRNGT (IINDA)
                  End If
               Else
!
!  Only B still with unprocessed values
!
                  IRNGT (IWRKD:IINDB) = JWRKT (IWRKD:IINDB)
                  IWRK = IWRKF
                  Exit
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         LMTNA = 2 * LMTNA
      End Do
!
!  Clean up
!
      Deallocate (JWRKT)
      Return
!
End Subroutine R_mrgref