C    order_bdm.f                May 2004  S.G. 
C
C  This program orders bdm catalogs by mass



      CHARACTER*120   name1,name2,f_tmp
      CHARACTER*120   line
      CHARACTER*90    trash
      PARAMETER (nmax = 120000)
      DIMENSION r_bdm(3,nmax),r_prim(3,3),
     +          v_fof(nmax),
     +          v_bdm3(3,nmax),v_bdm(nmax),
     +          total_mass_bdm(nmax),rh_bdm(nmax), Vrms_bdm(nmax),
     +          vcirc_bdm(nmax),npar_bdm(nmax),
     +          npar_bdm1(nmax),npar_bdm2(nmax),npar_bdm3(nmax),
     +          indx(nmax), rmax(nmax),conc(nmax)

      
c      box = 50.

      write(*,*) 'filename of the catshort '
      READ (*,'(A)') f_tmp         
c     f_tmp = 'Catshort.DAT'
      CALL  get_name(f_tmp,i,j)
      name2 = f_tmp(i:j)
      name1 = f_tmp(i:j)//'_ord'

      CALL  get_name(name1,in,jn)
      OPEN(11,FILE= name1(in:jn),FORM='FORMATTED',
     +       STATUS='UNKNOWN')

      CALL  get_name(name2,in,jn)
      OPEN(10,FILE= name2(in:jn),FORM='FORMATTED',
     +       STATUS='OLD')

      i = 1
      Do k = 1,21
        READ(10,'(A)')  trash
        write(11,'(A)')  trash
      ENDDO
      READ(10,'(A)')  line
      WRITE(11,123)  line
      write(*,123)  line
  123 Format('    #    ',A120)
      DO WHILE (.TRUE.)
      READ(10,*,END=260)
     +            r_bdm(1,i),r_bdm(2,i),r_bdm(3,i),
     +            v_bdm3(1,i),v_bdm3(2,i),v_bdm3(3,i),
     +            total_mass_bdm(i),rh_bdm(i), Vrms_bdm(i),
     +            vcirc_bdm(i),npar_bdm(i),
     +            rmax(i),conc(i)
      i = i + 1
      ENDDO
  260 n_bdm = i - 1
      write(*,*) n_bdm,'  objects read from the bdm_catalog'

      CALL  indexx(n_bdm,vcirc_bdm,indx)

      n_40 = 0
      n_50 = 0
      n_80 = 0
      n_100 = 0

      kk = 0
      DO k = n_bdm,1,-1
      kk = kk + 1
        i = indx(k)
        write(11,200) kk,
     +            r_bdm(1,i),r_bdm(2,i),r_bdm(3,i),
     +            v_bdm3(1,i),v_bdm3(2,i),v_bdm3(3,i),
     +            total_mass_bdm(i),rh_bdm(i), Vrms_bdm(i),
     +            vcirc_bdm(i),npar_bdm(i),
     +            rmax(i),conc(i)
      IF(vcirc_bdm(i).GT.100.) n_100 = n_100 + npar_bdm(i)
      IF(vcirc_bdm(i).GT.80.) n_80 = n_80 + npar_bdm(i)
      IF(vcirc_bdm(i).GT.50.) n_50 = n_50 + npar_bdm(i)
      IF(vcirc_bdm(i).GT.40.) n_40 = n_40 + npar_bdm(i)

      ENDDO

      write(*,*) n_100,' particles in halos with v_circ > 100 km/s'
      write(*,*)  n_80,' particles in halos with v_circ >  80 km/s'
      write(*,*)  n_50,' particles in halos with v_circ >  50 km/s'
      write(*,*)  n_40,' particles in halos with v_circ >  40 km/s'



 200  FORMAT(I8,3F8.3,3X,3F8.1,G12.3,3F8.2,I8,2F8.1)

      STOP
      END


******************************************************************* 
c     distance with periodic boundary conditions

      FUNCTION DIST(boxsize,XI,XJ,YI,YJ,ZI,ZJ)

      DX = ABS(XI - XJ)
      DDX =ABS(DX - boxsize)
      DX = MIN(DX,DDX)
      DY = ABS(YI - YJ)
      DDY =ABS(DY - boxsize)
      DY = MIN(DY,DDY)
      DZ = ABS(ZI - ZJ)
      DDZ =ABS(DZ - boxsize)
      DZ = MIN(DZ,DDZ)
      DIST = SQRT(DX**2 + DY**2 + DZ**2)
      RETURN
      END

*******************************************************************
C    Ordering of a real array
C
C    Indexing  Sect. 8.4. of Numerical Recipis
C    Changed to integer field

      SUBROUTINE indexx(n,arr,indx)
      INTEGER n,indx(n),M,NSTACK
      REAL arr(n)                       
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
      REAL a                            
      do 11 j=1,n
        indx(j)=j
11    continue
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 13 j=l+1,ir
          indxt=indx(j)
          a=arr(indxt)
          do 12 i=j-1,1,-1
            if(arr(indx(i)).le.a)goto 2
            indx(i+1)=indx(i)
12        continue
          i=0
2         indx(i+1)=indxt
13      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        itemp=indx(k)
        indx(k)=indx(l+1)
        indx(l+1)=itemp
        if(arr(indx(l+1)).gt.arr(indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l)).gt.arr(indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l+1)).gt.arr(indx(l)))then
          itemp=indx(l+1)
          indx(l+1)=indx(l)
          indx(l)=itemp
        endif
        i=l+1
        j=ir
        indxt=indx(l)
        a=arr(indxt)
3       continue
          i=i+1
        if(arr(indx(i)).lt.a)goto 3
4       continue
          j=j-1
        if(arr(indx(j)).gt.a)goto 4
        if(j.lt.i)goto 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
        goto 3
5       indx(l)=indx(j)
        indx(j)=indxt
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END
C  (C) Copr. 1986-92 Numerical Recipes Software .-35?421.1-9.

***********************************************************

      SUBROUTINE get_name(name,i,j)
      CHARACTER*120 name
      DO i=1,120
         IF(name(i:i).NE.' ') GOTO 1
      ENDDO
      i=1
      j=1
      RETURN
 1    CONTINUE
      DO j=i,120
         IF(name(j:j).EQ.' ') GOTO 2
      ENDDO
      j=121
 2    j=j-1
      IF(name(i:i).EQ.'''') i=i+1
      IF(name(j:j).EQ.'''') j=j-1
      END

