PROGRAM average6

c	This version uses direct access files rather than virtual memory because
c	of stupid system limitation on SG.

c	Average map

	PARAMETER	aa = 197.17, bb = 127.03, cc = 134.18, beta = 97.64
	PARAMETER	nx = 240, ny = 144, nz = 160		!2.8 Ang output
	PARAMETER	mx = 400, my = 256, mz = 270		!2.8 Ang input
	PARAMETER	nx2 = nx/2, ny2 = ny/2, nz2 = nz/2, nxm = nx-1
	PARAMETER	mx2 = mx/2, my2 = my/2, mz2 = mz/2, mxm = mx-1
	PARAMETER	nxl = 157  ! MAX(IFIX(nx*.65)+1,122)
	PARAMETER	nyl = 114  ! MAX(ny/2,114)
	PARAMETER	nzl = 120  ! MAX(IFIX(nz*.58)+1,120)

	EXTERNAL	time
	REAL		map(0:mxm,0:my2,0:mz2)
	INTEGER		in(5), maprec, averec, time
	CHARACTER*25	scrfil
	CHARACTER*80	title, filnam

	REAL		mat(3,3,5), fx(6), fy(6), fz(6)
	INTEGER*2	iout(3,6), ix, iy, iz

	REAL		aver(-nxl:nxl,-nzl:nzl)

	REAL*8		rms

c	Generate unique file name for scratch file based upon time in seconds

	scrfil(1:5) = 'AVER.'
	WRITE (scrfil(6:14), 19) time()
19	FORMAT (i9)

c	Read Ten Eyck map for asymmetric unit

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

	READ (1) title
	WRITE (6, 1001) title
1001	FORMAT (1x,a)

	DO iy = 0, my2
	  READ (1) in
c	  WRITE (6, 1002) in
1002	  FORMAT (5i5)
	  READ (1) ((map(ix, iy, iz), iz=0,mz2), ix= 0,mxm)
	endDO

	READ (1) in
c	  WRITE (6, 1002) in

	CLOSE (unit=1)

c	Read local symmetry operations

	WRITE (6, 1000) 'Enter file name for local symmetry operations:'
	READ (5, 1000) filnam
	OPEN (unit=1, status='OLD', form='FORMATTED', readonly, file=filnam)

	READ (1, 1003)
	DO n = 1,5
	  READ (1, 1003) ((mat(i,j,n), i=1,3), j=1,3)
1003	  FORMAT (9f10.7)
	endDO

	CLOSE (unit=1)

c	Open file for list of equivalent points

	WRITE (6, 1000) 'Enter file name for list of equivalent points:'
	READ (5, 1000) filnam
	OPEN (unit=1, status='OLD', form='UNFORMATTED', readonly, file=filnam)

c	Loop through point list

c	NOTE:  Integers are coordinates in asymmetric unit while reals
c	       are coordinates in region where local symmetry applies

	nin = 0
	averec = -nyl
	DO ii = -nxl, nxl
	  DO jj = -nzl, nzl
	    aver(ii,jj) = 0.0
	  endDO
	endDO

	OPEN (unit=11, access='DIRECT', recl=(2*nxl+1)*(2*nzl+1),
	1	status='NEW', file=scrfil, maxrec=2*nyl+1)

	DO WHILE (.TRUE.)

	  READ (1, end=9) ix, iy, iz, fx(1), fy(1), fz(1)	  
	  nin = nin + 1
	  IF (MOD(nin,10000) .EQ. 0) THEN
c	    WRITE (6, 5040) nin	
5040	    FORMAT (' Averaging point', i10)
	  endIF

c	Apply local symmetry

	  DO n = 1,5
	    fx(n+1) = fx(1)*mat(1,1,n) + fy(1)*mat(1,2,n) + fz(1)*mat(1,3,n)
	    fy(n+1) = fx(1)*mat(2,1,n) + fy(1)*mat(2,2,n) + fz(1)*mat(2,3,n)
	    fz(n+1) = fx(1)*mat(3,1,n) + fy(1)*mat(3,2,n) + fz(1)*mat(3,3,n)
	  endDO

c	Since points are sorted in increasing sections of y, if iy > averec
c	write out section and work on next one

	  DO WHILE (iy. GT. averec) 
	    WRITE (11, rec=averec+nyl+1) aver
	    DO ii = -nxl, nxl
	      DO jj = -nzl, nzl
	        aver(ii,jj) = 0.0
	      endDO
	    endDO
	    averec = averec+1
	  endDO

c	Put in asymmetric unit

	  DO n = 1,6
	    fy(n) = MOD(fy(n)+10.,1.)

	    IF (fy(n) .GT. 0.5) THEN
	      fx(n) = MOD(fx(n)+10.5,1.)
	      fz(n) = MOD(fz(n)+10.5,1.)
	      fy(n) = fy(n) - 0.5
	    ELSE
	      fx(n) = MOD(fx(n)+10.,1.)
	      fz(n) = MOD(fz(n)+10.,1.)
	    endIF

	    IF (fz(n) .GT. 0.5) THEN
	      fz(n) = 1.0 - fz(n)
	      IF (fx(n) .NE. 0.0) fx(n) = 1.0 - fx(n)
	    endIF

	    ixp = MOD(NINT(fx(n)*mx),mx)
	    iyp = MOD(NINT(fy(n)*my),my)
	    izp = MOD(NINT(fz(n)*mz),mz)

	    aver(ix,iz) = aver(ix,iz) + map(ixp,iyp,izp)
	  endDO

	  aver(ix,iz) = aver(ix,iz)/6.

	endDO

9	CLOSE (unit=1)

	DO WHILE (nyl .GE. averec) 
	  WRITE (11, rec=averec+nyl+1) aver
	  DO ii = -nxl, nxl
	    DO jj = -nzl, nzl
	      aver(ii,jj) = 0.0
	    endDO
	  endDO
	  averec = averec+1
	endDO

c	Calculate rms value of map and apply dyad

	npts = 0
	rms = 0.0d0
	DO iy = -nyl, nyl
	  READ (11, rec=iy+nyl+1) aver
	  DO iz = 0, nzl
	    DO ix = -nxl, nxl
	      aver(-ix,-iz) = aver(ix,iz)
	      IF (aver(ix,iz) .NE. 0.0) THEN
	        rms = rms + aver(ix,iz)*aver(ix,iz)
	        npts = npts + 1
	      endIF
	    endDO
	  endDO
	  WRITE (11, rec=iy+nyl+1) aver
	endDO

	rms = SQRT(rms/npts)

	WRITE (6, 1004) npts, rms
1004	FORMAT (' Rms of', i7, ' points in averaged map:', f10.4)
	
c	Write averaged map in 2 halves

	WRITE (6, 1000) 'Enter title for map:'
	READ (5, 1000) title
	WRITE (6, 1000) title

	WRITE (6, 1000) 'Enter file name for first half of averaged map:'
	READ (5, 1000) filnam
	WRITE (6, 1000) filnam
	OPEN (unit=1, status='UNKNOWN', form='UNFORMATTED', file=filnam)

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

	in(2) = nz-nzl
	in(3) = nz+nzl
	in(4) = nx-nxl
	in(5) = nx+nxl

	DO iy = -nyl, 0
	  READ (11, rec=iy+nyl+1) aver
	  in(1) = ny+iy
	  WRITE (1) in
c	  WRITE (6, 1002) in
	  WRITE (1) ((aver(ix,iz), iz=-nzl,nzl), ix=-nxl,nxl)
	endDO

	WRITE (1) -1, -1, -1, -1, -1

	CLOSE (unit=1)

	WRITE (6, 1000) 'Enter file name for second half of averaged map:'
	READ (5, 1000) filnam 
	WRITE (6, 1000) filnam
	OPEN (unit=1, status='UNKNOWN', form='UNFORMATTED', file=filnam) 

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

	in(2) = nz-nzl
	in(3) = nz+nzl
	in(4) = nx-nxl
	in(5) = nx+nxl

	DO iy = 0, nyl
	  READ (11, rec=iy+nyl+1) aver 
	  in(1) = ny+iy 
	  WRITE (1) in 
c	  WRITE (6, 1002) in
	  WRITE (1) ((aver(ix,iz), iz=-nzl,nzl), ix=-nxl,nxl)
	endDO

	WRITE (1) -1, -1, -1, -1, -1

	CLOSE (unit=1)
	CLOSE (unit=11, disp='DELETE')

	end