        program historgram

          implicit none

          real*4 r_data(32768)
          integer i_cdata(1000)
          integer i_rdata(1000)
          integer i_minrc(100000)
          integer i_maxrc(100000)
          integer i_minr
          integer i_maxr
          integer i_minc
          integer i_maxc
          integer i_pts

          real*4 r_hist(0:256)
          real*4 r_nullv
          real*4 r_dist
          real*4 r_minv
          real*4 r_maxv


          character*255 a_datafile
          character*255 a_pntsfile
	  character*255 a_value
          character*512 a_cmnd
	  
          integer i
          integer j
          integer k
          integer i_c
          integer i_r
          integer i_h

	  integer i_arg
          integer i_inarg
          integer i_samps
          integer i_sfmt
          integer i_shdr
          integer i_rhdr
          integer i_chdr
  
          integer iargc
          external iargc

          integer length
          external length
  
          i_sfmt = 4
          
c          write(6,*) 'initializing min/max rc'
          do i=1,100000
            i_minrc(i)=-1
            i_maxrc(i)=-1
          end do

          do i=0,256
            r_hist(i)=0.0
          end do

c          write(6,*) 'getting command line arguments'
          i_inarg = iargc()
          if (i_inarg .lt. 4) then
            call write_greeting()
            stop 'done'
          else
            call getarg(1,a_pntsfile)
            call getarg(2,a_datafile)
            call getarg(3,a_value)
            read(a_value,*) i_sfmt
            write(6,*) 'i_setvfmt = ',i_sfmt
            call getarg(4,a_value)
            read(a_value,*) i_samps
            write(6,*) 'i_setcols = ',i_samps
            if (i_inarg .ge. 5) then
              call getarg(5,a_value)
              read(a_value,*) i_shdr
              write(6,*) 'i_setshdr = ',i_shdr
            else
              i_shdr = 0
            end if
            if (i_inarg .ge. 6) then
              call getarg(6,a_value)
              read(a_value,*) i_rhdr
              write(6,*) 'i_setrhdr = ',i_rhdr
            else
              i_rhdr = 0
            end if
            if (i_inarg .ge. 7) then
              call getarg(7,a_value)
              read(a_value,*) i_chdr
              write(6,*) 'i_setchdr = ',i_chdr
            else
              i_rhdr = 0
            end if
            if (i_inarg .ge. 8) then
              call getarg(8,a_value)
              read(a_value,*) r_nullv
            else
              r_nullv=-1.e20
            end if
          end if

          if (i_sfmt .ne. 4) stop '*** Sample format not supported ***'
          
          write(6,*) 'record length in bytes = ',i_samps*4+i_rhdr
          open(unit=20,file=a_datafile,status='old',form='unformatted',access='direct',recl=i_samps*4+i_rhdr)
          open(unit=30,file=a_pntsfile,status='old',form='formatted')
          open(unit=40,file='histogram.dat',status='unknown',form='formatted')

          i_pts=0

c          write(6,*) 'reading in points'          
          i_minr=1e10
          i_maxr=-1e10
          i_minc=1e10
          i_maxc=-1e10
          do i=1,1000
            read(30,*,err=900,end=900) i_cdata(i),i_rdata(i)
            i_pts=i
            i_minr=min(i_minr,i_rdata(i))
            i_maxr=max(i_maxr,i_rdata(i))
            i_minc=min(i_minc,i_cdata(i))
            i_maxc=max(i_maxc,i_cdata(i))
        
          end do
          write(6,*) 'Too many points - only using first 1000'
900       continue

          write(6,*) 'min/max row = ',i_minr,i_maxr
          write(6,*) 'min/max col = ',i_minc,i_maxc

          if (i_maxr-i_minr .ge. 1000000) stop 'Region too big. Must be less than 1000000 rows'

          do i=1,i_pts
c            write(6,*) 'at point: ',i
            j=i+1
            if (j .gt. i_pts) j=1
            r_dist = sqrt(float(i_cdata(j)-i_cdata(i))**2.+float(i_rdata(j)-i_rdata(i))**2.)
c            write(6,*) ' ',i_rdata(i),i_cdata(i),r_dist
            do k=0,r_dist*4
              i_c = i_cdata(i)+(i_cdata(j)-i_cdata(i))*k/(4*r_dist)
              i_r = i_rdata(i)+(i_rdata(j)-i_rdata(i))*k/(4*r_dist)
c              write(6,*) '   ',i_r,i_c
              if (i_minrc(i_r-i_minr+1) .eq. -1) i_minrc(i_r-i_minr+1)=i_c
              if (i_maxrc(i_r-i_minr+1) .eq. -1) i_maxrc(i_r-i_minr+1)=i_c
              i_minrc(i_r-i_minr+1) = min(i_minrc(i_r-i_minr+1),i_c)
              i_maxrc(i_r-i_minr+1) = max(i_maxrc(i_r-i_minr+1),i_c)
            end do
          end do

          r_minv=+1e27
          r_maxv=-1e27
          do i_r=0,i_maxr-i_minr
c            write(6,*) 'Reading line: ',i_r+i_minr+1,'   Samps: ',i_minrc(i_r+1)+1,i_maxrc(i_r+1)+1
            read(20,rec=i_r+i_minr+1+i_shdr/(i_samps*4+i_rhdr)) (r_data(i_c),i_c=1,i_rhdr/4+min(i_samps,i_maxrc(i_r+1)+1))
            do i_c=i_minrc(i_r+1)+1,i_maxrc(i_r+1)+1
              if (r_data(i_c+i_rhdr/4) .gt. r_nullv) then
                r_minv=min(r_minv,r_data(i_c+i_rhdr/4))
                r_maxv=max(r_maxv,r_data(i_c+i_rhdr/4))
              end if
            end do
          end do

          write(6,*) 'min/max data values: ',r_minv,r_maxv

          i_pts = 0
          do i_r=0,i_maxr-i_minr
            read(20,rec=i_r+i_minr+1+i_shdr/(i_samps*4+i_rhdr)) (r_data(i_c),i_c=1,i_rhdr/4+min(i_samps,i_maxrc(i_r+1)+1))
            do i_c=i_minrc(i_r+1)+1,i_maxrc(i_r+1)+1
              if (i_pts .lt. 999999) then
                r_hist(nint(256*(r_data(i_c+i_rhdr/4)-r_minv)/(r_maxv-r_minv))) = r_hist(nint(256*(r_data(i_c+i_rhdr/4)-r_minv)/(r_maxv-r_minv))) + 1
                i_pts = i_pts +1
              end if
            end do
          end do

          do i_h = 0,256
            r_hist(i_h) = r_hist(i_h)/i_pts
            write(40,*) i_h*(r_maxv-r_minv)/256.+r_minv,r_hist(i_h)
          end do
          write(a_cmnd,'(a,i6,a)') 'xmgrace -geometry 800x600 -noask -free -pexec "s0 line type 3" -pexec "s0 dropline on" '//
     &         '-pexec ''yaxis label "Fraction of Total"''  -pexec ''xaxis label "Value"''  -pexec ''title "'//
     &         a_datafile(1:length(a_datafile))//'"''  '//'-pexec ''subtitle "Total number of points:',i_pts,'"''  '//
     &         ' histogram.dat &'
          call system(a_cmnd)
          write(6,*) 'Histogram Done'
        end

c****************************************************************

      integer*4 function length(a_string)

c****************************************************************
c**   
c**   FILE NAME: rdf_reader.f
c**   
c**   DATE WRITTEN: 15-Sept-1997
c**   
c**   PROGRAMMER: Scott Shaffer
c**   
c**   FUNCTIONAL DESCRIPTION: This function returns the position 
c**   of the last none blank character in the string. 
c**   
c**   ROUTINES CALLED:
c**   
c**   NOTES: 
c**   
c**   UPDATE LOG:
c**   
c**   Date Changed        Reason Changed                  CR # and Version #
c**   ------------       ----------------                 -----------------
c**   
c*****************************************************************

      implicit none

c     INPUT VARIABLES:

      character*(*) a_string
      
c     OUTPUT VARIABLES:

c     LOCAL VARIABLES:

      integer i_len

c     COMMON BLOCKS:

c     EQUIVALENCE STATEMENTS:

c     DATA STATEMENTS:

c     FUNCTION_STATEMENTS:

c     PROCESSING STEPS:

c      write(6,*) 'here =',a_string(1:60)
      i_len=len(a_string)
      do while(i_len .gt. 0 .and. (a_string(i_len:i_len) .eq. ' ' .or. 
     &     ichar(a_string(i_len:i_len)) .eq. 0))
         i_len=i_len-1
c         write(6,*) i_len,' ',ichar(a_string(i_len:i_len)),' ',a_string(i_len:i_len)
      enddo

      length=i_len
      return

      end


        subroutine write_greeting()
          implicit none

            write(6,*) 'This is a very unfriendly program.  Figure it out for youself'

          return
        end
