PROGRAM corrtemap

c	Calculates correlation between two maps.  Program ignores those map
c	points which are zero.

	REAL*4		title(20)
	REAL		map1(200000), map2(200000)
	INTEGER		in(5)
	CHARACTER*80	filnam

	WRITE (6, 1000) 'Enter file name of first Ten Eyck map:'
1000	FORMAT (a)
	READ (5, 1000) filnam
	OPEN (unit=1, file=filnam, form='UNFORMATTED', status='OLD', readonly)

	WRITE (6, 1000) 'Enter file name of second Ten Eyck map:'
	READ (5, 1000) filnam
	OPEN (unit=2, file=filnam, form='UNFORMATTED', status='OLD', readonly)

c	Read titles

	READ (1) title
	WRITE (6, 1001) title
1001	FORMAT (20a4)
	READ (2) title
	WRITE (6, 1001) title

c	Loop through sections

	DO WHILE (.TRUE.)
	  READ (1) in
	  READ (2) in

	  IF (in(1) .LT. 0) GOTO 9

	  npts = (in(5)-in(4)+1)*(in(3)-in(2)+1)
	  READ (1) (map1(i), i=1,npts)
	  READ (2) (map2(i), i=1,npts)

	  sum1 = 0.d0
	  sum2 = 0.d0
	  nonzero = 0
	  DO i = 1,npts
	    IF (map1(i)*map2(i) .NE. 0.0) THEN
	      sum1 =  map1(i) + sum1
	      sum2 =  map2(i) + sum2
	      nonzero = nonzero + 1
	    endIF
	  endDO
	  av1 = sum1/nonzero
	  av2 = sum2/nonzero
	  sum1 = 0.d0
	  sum2 = 0.d0
	  sum3 = 0.d0
	  DO i = 1,npts
	    IF (map1(i)*map2(i) .NE. 0.0) THEN
	      map1(i) = map1(i)-av1
	      map2(i) = map2(i)-av2
	      sum1 = sum1 + map1(i)*map1(i)
	      sum2 = sum2 + map2(i)*map2(i)
	      sum3 = sum3 + map1(i)*map2(i)
	    endIF
	  endDO
	  corr = sum3/(SQRT(sum1*sum2))
	  WRITE (6, 1002) in, av1, av2, corr, nonzero
1002	  FORMAT (5i5, 3f15.4, i10)
	endDO

9	CLOSE (unit=1)
	CLOSE (unit=2)

	end