program statmat
      character*70 filename(5),statfile
      character*5 done
      integer*4 col(5),row(5),count
      real*4 nobss,data(500,500,5),colat(5),long(5),space(5),
     >       xmean,ymean
      dimension xdata(250000),ydata(250000)
c
c-------------------------------- program description
c	statmat defines the basic statistics of the input grids.  see
c	the write statements for the specific values calculated.  also,
c	the code loops through all input grids and calculates the
c	correlation coefficients between all combinations of input
c	data.
c
c
      write (*,*) 'OUTPUT STATISTICS FILE'
      read (*,9990) statfile
      open (25, file=statfile, form='formatted')
c
c----------------read all data into array (row,col,layer or map number)
      i=1
 10   write (*,*) 'INPUT MATRIX WITH TPLOT HEADER'
      read (*,9990) filename(i)
 9990 format (a70)
      open (10, file=filename(i),status='old',form='formatted')
c
      read (10,*) col(i)
      read (10,*) row(i)
      read (10,*) colat(i)
      read (10,*) long(i)
      read (10,*) space(i)
      read (10,*) ((data(j,k,i),k=1,col(i)),j=1,row(i))
      close (10)
      write (*,*) 'ARE YOU THROUGH YET???'
      read (*,9991) done
 9991 format (a5)
      if (done .eq. 'y'  .or.  done .eq. 'yes') go to 50
      i=i+1
      go to 10
c------------------------maximum number of input data sets =count
c                        loops from line 80 to 400 increment through all 
c                        maps comparing all possible logical 
c                        combinations of map to map
 50   count=i
      write (*,*) 'total number of input data sets =',count
      do 60 i=1,count-1
         if (col(i) .ne. col(i+1)  .or.  row(i) .ne. row(i+1)) then
            write (*,*) filename(i),col(i),row(i)
            write (*,*) filename(i),col(i),row(i)
            write (*,*) 'rows or columns do not match'
            stop 0001
         endif
 60   continue
c
c-------------------------from 200 to write statement of variables is
c                         the statistical calculations using two
c                         references:
c                          1) Davis, Statistics and Data Analysis in 
c                             Geology, 2nd ed., 1986 pp. 41
c                          2) Young, Statistical Treatment of Experi-
c                             mental Data, 1962, McGraw Hill, 115-132
c
 80   continue
      do 400 icnt=1,count
         do 400 jcnt=icnt,count
c
            do 100 j=1,row(icnt)
               do 100 k=1,col(icnt)
                  ij=(col(icnt)*(j-1))+k
                  xdata(ij)=data(j,k,icnt)
                  ydata(ij)=data(j,k,jcnt)
 100        continue 
c-------------------------loops that sum x, x**2, y, y**2 and xy
      nobs=row(icnt)*col(icnt)
      if (nobs .ne. ij) stop 0002
      nobss=float(nobs)
      xsum=0.0
      xsumsqr=0.0
      ysum=0.0
      ysumsqr=0.0
      sumxy=0.0
      xmax=xdata(1)
      xmin=xmax
      ymax=ydata(1)
      ymin=ymax
      do 240 j4=1,nobs
         xsum=xsum+xdata(j4)
         xsumsqr=xsumsqr+(xdata(j4))**2
         ysum=ysum+ydata(j4)
         ysumsqr=ysumsqr+(ydata(j4))**2
         sumxy=sumxy+(xdata(j4)*ydata(j4))
         xmax=max(xmax,xdata(j4))
         xmin=min(xmin,xdata(j4))
         ymax=max(ymax,ydata(j4))
         ymin=min(ymin,ydata(j4))
 240  continue
c------------------------find corrected sum of products, covariance
c                        and corrected sum of squares (x) (y)
c
      xmean=xsum/nobss
      ymean=ysum/nobss
      sumprod=sumxy-((xsum*ysum)/nobss)
      covarxy=sumprod/(nobss-1.0)
      xcsumsqr=xsumsqr-((xsum**2)/nobss)
      ycsumsqr=ysumsqr-((ysum**2)/nobss)
c
c------------------------find variance, standard deviation for x and y
c     
      xvar=xcsumsqr/(nobss-1.0)
      yvar=ycsumsqr/(nobss-1.0)
      xstdev=sqrt(xvar)
      ystdev=sqrt(yvar)
c-----------------------find correlation coefficient by Davis method
      corrDxy=covarxy/(xstdev*ystdev)
c-----------------------find slopes, intercepts and correlation 
c                       coefficient by Young method
      xslope=((nobss*sumxy)-(xsum*ysum))/((nobss*xsumsqr)-xsum**2)
      yslope=((nobss*sumxy)-(xsum*ysum))/((nobss*ysumsqr)-ysum**2)
      xintcpt=((ysum*xsumsqr)-(sumxy*xsum))/((nobss*xsumsqr)-xsum**2)          
      yintcpt=((xsum*ysumsqr)-(sumxy*ysum))/((nobss*ysumsqr)-ysum**2)
      corrYxy=sqrt(xslope*yslope)
c
c------------------------write out this mess for individual pass and 
c                        overlapping lengths of passes
c
      write (25,*) 'X = ',filename(icnt)
      write (25,*) 'Y = ',filename(jcnt)
      write (25,9992) xmean,ymean,xvar,yvar,xstdev,ystdev
 9992 format('X MEAN =',f9.3,' Y MEAN =',f9.3,/,
     >       'X VARIANCE=',f9.1,' Y VARIANCE=',f9.1,' X STDEV=',
     >        f8.3,' Y STDEV=',f8.3)
      write (25,9993) covarxy,corrDxy
 9993 format ('COVARIANCE XY=',f9.1,' Davis CORRELATION COEF=',f8.3)
      write (25,9994) xslope,xintcpt,yslope,yintcpt,corrYxy
 9994 format ('X SLOPE=',f8.3,' X INTERCEPT=',f8.3,' Y SLOPE=',
     >         f8.3,' Y INTERCEPT=',f8.3,/,'Young CORRELATION COEF=',
     >         f8.3,)
      write (25,9995) xmax,xmin,ymax,ymin
      write (*,*) corrDxy
 9995 format('X-MAX=',f9.3,' X-MIN=',f9.3,' Y-MAX=',f9.3,
     >       ' Y-MIN=',f9.3,/)
c
c
c-----------------------increment to next set of passes
 400  continue
c
 999  continue
      close (10)
      close (25)
      stop
      end