PROGRAM		corel

c	Calculate corelation between 2 maps

	PARAMETER	nx = 240, ny = 144, nz = 160		!2.8 Ang output
	PARAMETER	nx1 = 0,   nx2 = 100
	PARAMETER	ny1 = -30, ny2 = 40
	PARAMETER	nz1 = -10, nz2 = 85

	INTEGER		in(5)
	CHARACTER*80	title
	CHARACTER*40	filename

	REAL		map(2,nx1:nx2,ny1:ny2,nz1:nz2)
	REAL*8		sum1, sum2, sum11, sum12, sum22

c	Read maps

	DO nmap = 1,2
1	  WRITE (6, 1001) '$Enter file name for map', nmap, ':     '
1001	  FORMAT (a, i3, a)
	  READ (5, 1001, end=10) filename
	
	  OPEN (unit=1, status='OLD', form='UNFORMATTED',
	1	file=filename, err=1)

	  READ (1) title
c	  WRITE (6, 1001) title

	  DO WHILE (.TRUE.)
	    READ (1) in
c	    WRITE (6, 1002) in
1002	    FORMAT (5i5)
	    IF (in(1) .EQ. -1) GOTO 5

	    in(1) = in(1)-nz
	    in(2) = in(2)
	    in(3) = in(3)
	    in(4) = in(4)-ny
	    in(5) = in(5)-ny

	    READ (1) ((map(nmap,ix,iy,in(1)), ix=in(2),in(3)), iy=in(4),in(5))

	  endDO

5	  CLOSE (unit=1)
	endDO

c	Calculate corellation coefficient

10	DO i = nx1,nx2
	  DO j = ny1,ny2
	    DO k = nz1,nz2
	      sum1 = sum1 + map(1,i,j,k)
	      sum2 = sum2 + map(2,i,j,k)
	      sum11 = sum11 + map(1,i,j,k)*map(1,i,j,k)
	      sum12 = sum12 + map(1,i,j,k)*map(2,i,j,k)
	      sum22 = sum22 + map(2,i,j,k)*map(2,i,j,k)
	   endDO
	  endDO
	endDO

	npts = (nx2-nx1+1)*(ny2-ny1+1)*(nz2-nz1+1)
	sum1 = sum1/npts
	sum2 = sum2/npts
	sum11 = sum11/npts
	sum12 = sum12/npts
	sum22 = sum22/npts
	sig1 = SQRT(sum11-sum1*sum1)
	sig2 = SQRT(sum22-sum2*sum2)
	corr = (sum12-sum1*sum2)/(sig1*sig2)

	WRITE (6, 1001)
	WRITE (6, 1004) sum1, sig1, sum2, sig2, corr
1004	FORMAT (' Mean 1 =', t10, f10.3, t30, ' Sigma 1 =', t40, f10.3, /,
	1	' Mean 2 =', t10, f10.3, t30, ' Sigma 2 =', t40, f10.3, /,
	2	' Corr =', t10, f10.5)

	end