cProgram to perform some sanity checks on the initial conditiona
      implicit double precision (a-h,o-z)
      real msun
      parameter (msun = 1.989e33, pc = 3.085678e18, 
     &     nmax=30000,ns=12288, g = 6.670e-8*msun*1e-10/pc)
      double precision mass(nmax),totmass,temp(nmax),temp2(nmax),
     &     temp3(nmax),temp4(nmax)
      dimension r(nmax,3),v(nmax,3),rbar(3),vbar(3),rrel(3),vrel(3)
      real m(2),h(3)
      data totmass/0.d0/,rbar/3*0.d0/,vbar/3*0.d0/
      open (7,file='ic.dat')
      i = 1
 10   continue
      read (7,*,end=20) mass(i),(r(i,j),j=1,3),(v(i,j),j=1,3)
      totmass = totmass + mass(i)
      do j = 1,3
         rbar(j) = rbar(j) + mass(i)*r(i,j)
         vbar(j) = vbar(j) + mass(i)*v(i,j)
      enddo
      i = i + 1
      goto 10
 20   continue
      n = i-1
      write (6,*) 'total mass, N ',totmass, n
      write (6,*) 'cm position ',(rbar(j)/totmass,j=1,3)
      write (6,*) 'cm velocity ',(vbar(j)/totmass,j=1,3)
      nb = (n - ns)/2
cCheck median masses
      do i = 1,ns
         temp(i) = mass(i)
      enddo
      x1 = select(ns/2,ns,temp)
      x2 = select(ns/2+1,ns,temp)
      write (6,*) 'median single mass ', 0.5*(x1 + x2)
      do i = 1,nb
         temp(i) = mass(ns + 2*i - 1) + mass(ns + 2*i)
      enddo
      x1 = select(nb/2,nb,temp)
      x2 = select(nb/2+1,nb,temp)
      write (6,*) 'median binary mass ', 0.5*(x1 + x2)
cCheck median radius
      do i = 1,n
         temp(i) = r(i,1)**2 + r(i,2)**2 + r(i,3)**2 
      enddo
      x1 = select(n/2,n,temp)
      x2 = select(n/2+1,n,temp)
      write (6,*) 'median radius ', 0.5*(sqrt(x1) + sqrt(x2))
cTest binaries
c1. Distribution of q
      do i = 1,nb
         temp(i) = mass(ns+2*i-1)/mass(ns+2*i)
      enddo
      call sort(nb,temp)
      open(9,file='sorted_q')
      do i = 1,nb
         write (9,*) i,temp(i)
      enddo
c2. Distribution of log a, e2
      do i = 1,nb
         m(1) = mass(ns+2*i-1)
         m(2) = mass(ns+2*i)
         r2 = 0.0
         v2 = 0.0
         rdotv = 0.0
         do k = 1,3
            rrel(k) = r(ns+2*i-1,k) - r(ns+2*i,k)
            vrel(k) = v(ns+2*i-1,k) - v(ns+2*i,k)
            r2 = r2 + rrel(k)**2
            v2 = v2 + vrel(k)**2
         enddo
         h(1) = rrel(2)*vrel(3) - rrel(3)*vrel(2)
         h(2) = rrel(3)*vrel(1) - rrel(1)*vrel(3)
         h(3) = rrel(1)*vrel(2) - rrel(2)*vrel(1)
         energy = 0.5*v2 - g*(m(1)+m(2))/sqrt(r2)
         a = - g*(m(1)+m(2))/(2*energy)
         h2 = h(1)**2 + h(2)**2 + h(3)**2
         e2 = 1 - h2/(g*(m(1)+m(2))*a)
         temp(i) = log10(a)
         temp2(i) = e2
         temp3(i) = h(3)/sqrt(h2)
         temp4(i) = atan2(h(1),h(2))
      enddo
      call sort(nb,temp)
      call sort(nb,temp2)
      call sort(nb,temp3)
      call sort(nb,temp4)
      write (6,*) 'median a, e ', 10**temp(nb/2), sqrt(temp2(nb/2))
      open(10,file='sorted_loga')
      open(11,file='sorted_e2')
      open(12,file='sorted_cosi')
      open(13,file='sorted_node')
      do i = 1,nb
         write (10,*) i,temp(i)
         write (11,*) i,temp2(i)
         write (12,*) i,temp3(i)
         write (13,*) i,temp4(i)
      enddo
      stop
      end
c
c
      SUBROUTINE SORT(N,RA)
      double precision ra
      DIMENSION RA(N)
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          RRA=RA(L)
        ELSE
          RRA=RA(IR)
          RA(IR)=RA(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            RA(1)=RRA
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(RA(J).LT.RA(J+1))J=J+1
          ENDIF
          IF(RRA.LT.RA(J))THEN
            RA(I)=RA(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        RA(I)=RRA
      GO TO 10
      END
c
c
      FUNCTION select(k,n,arr)
      INTEGER k,n
      double precision select,arr(n)
      INTEGER i,ir,j,l,mid
      REAL a,temp
      l=1
      ir=n
1     if(ir-l.le.1)then
        if(ir-l.eq.1)then
          if(arr(ir).lt.arr(l))then
            temp=arr(l)
            arr(l)=arr(ir)
            arr(ir)=temp
          endif
        endif
        select=arr(k)
        return
      else
        mid=(l+ir)/2
        temp=arr(mid)
        arr(mid)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        if(j.ge.k)ir=j-1
        if(j.le.k)l=i
      endif
      goto 1
      END
C  (C) Copr. 1986-92 Numerical Recipes Software ]vrt43D04.