      implicit NONE
      
      integer max
      parameter (max=120000000)

      integer j,ndat
      real*8 diam(max),buf
      real*4 xmin,xmax,ymin,ymax,xx(max),yy(max),xval,yval

c     read data
      open(1,file='test_sorted',status='old')
      j=1
 100  continue
      read(1,*,end=200)buf,buf,buf,buf,diam(j)
      j=j+1
      if(j.gt.max) goto 200
      goto 100
 200  continue
      close(1)
      ndat=j-1
      write(*,*)'ndat=',ndat
      
      call pgbeg(0,'?',1,1) 
      call pgpap(0.,1.5)
      call pgsch(2.0)
      call pgslw(4)
    
      xmin=log10(0.01)
      xmax=log10(10.0)
      ymin=0.0
      ymax=8.2

      call pgvport(0.2,0.95,0.2,0.9)
      call pgswin(xmin, xmax, ymin, ymax)
      call pgbox('blcnst', 0., 0, 'bclnst', 0., 0)
      call pglab('\frDiameter D (km)', 
     &     '\frNumber of NEOs larger than D', ' ')  

      call pgsls(4)
      call pgslw(1)
      do j=1,21
         xval=-4+j
         call pgmove(xval,ymin)
         call pgdraw(xval,ymax)
      end do
      do j=1,21
         xval=log10(0.0003)+j
         call pgmove(xval,ymin)
         call pgdraw(xval,ymax)
      end do
      do j=1,9
         yval=0+j
         call pgmove(xmin,yval)
         call pgdraw(xmax,yval)
      end do
      do j=1,9
         yval=log10(0.3)+j
         call pgmove(xmin,yval)
         call pgdraw(xmax,yval)
      end do
      call pgsls(1)

      call pgslw(4)
c      call shell(ndat,hmag)
      do j=1,ndat
         xx(j)=log10(diam(j))
         yy(j)=log10(float(j))
      end do
      call pgline(ndat,xx,yy)

      call pgend


      end

c---------------------------------------------------------------------

      subroutine shell(n,a)
      integer n
      real*8 a(n)
      integer i,j,inc
      real*8 v
      inc=1 
 1    inc=3*inc+1
      if(inc.le.n) goto 1
 2    continue
      inc = inc/3
      do i=inc+1,n
         v = a(i)
         j = i
 3       if(a(j-inc).gt.v)then
            a(j)=a(j-inc)
            j=j-inc
            if(j.le.inc) goto 4
            goto 3
         end if
 4       a(j)=v
      end do
      if(inc.gt.1) goto 2
      return
      END
