PROGRAM distr
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      parameter (ndp0=1200)
      DIMENSION probd(ndp0)
      character*80 fin,fout,st
C     
      WRITE(6,*)'distribution analysis'
      iper=0
      call getarg(1,fin)
      call getarg(2,fout)
      if(iargc().gt.2)then
       call getarg(3,st)
       read(st,*)nag
       call getarg(4,st)
       read(st,*)ndp
       call getarg(5,st)
       read(st,*)xmin
       call getarg(6,st)
       read(st,*)xmax
       if(iargc().gt.6)then
        call getarg(7,st)
        read(st,*)per
        iper=1
        write(6,*)'periodic correction'
       endif
      else
       nag=1
       ndp=100
       xmin=-180.0d0
       xmax= 180.0d0
      endif

      dd=(xmax-xmin)/dble(ndp)

      if(ndp.gt.ndp0)call report('too many points')
      do 10 i=1,ndp
10    probd(i)=0.0d0
      an=0.0d0

      open(9,file=fin)
1     read(9,*,end=999,err=999)(par,i=1,nag)
      an=an+1.0d0
      if(par.le.xmax.and.par.ge.xmin)call digest(par,xmin,dd,probd,ndp)
      if(iper.eq.1)then
c      try to add and subtract the period, to have smooth ends
       pp=par+per
       if(pp.le.xmax.and.pp.ge.xmin)call digest(pp,xmin,dd,probd,ndp)
       pp=par-per
       if(pp.le.xmax.and.pp.ge.xmin)call digest(pp,xmin,dd,probd,ndp)
      endif
      goto 1
999   close(9)

      open(33,file=fout)
      dd=(xmax-xmin)/dble(ndp)
      fac=1.0d0/an/dd
      do 17 j=1,ndp
17    write(33,3333)xmin+dd*dble(2*j-1)/2.0d0,probd(j)*fac
3333  format(2e12.4)
      close(33)
      
      write(6,6009)int(an)
6009  format(i7,' points')

      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine digest(par,xmin,dd,probd,ndp)
      implicit none
      real*8 par,xmin,dd,probd(*),c,v,wp,w
      integer*4 j,jp,ndp
      j=int((par-xmin)/dd)+1
c     center of the grid cell:
      c=xmin+(dble(j)-0.5d0)*dd
      v=par-c
      if(v.gt.0)then
       jp=j+1
       wp=v/dd
      else
       jp=j-1
       wp=-v/dd
      endif
      w=1.0d0-wp
      if( j.ge.1.and.j .le.ndp)probd(j )=probd(j )+w
      if(jp.ge.1.and.jp.le.ndp)probd(jp)=probd(jp)+wp
      return
      end