C*VARIMAP.FOR******************************************************************** C * C This program calculates MEAN AND SD ARRAYS map from a set of density C maps. C The variance is calculated from a running average of the input * C * C Two maps of the same dimensions as the input map are output * C * C (i) Standard deviation = SQRT(variance(x,y)) * C (ii) real space average * C C Input maps must be calculated on the same grid. The user is * C responsible for scaling between maps. Input cards consist of * C a title card which is saved as the first title of the map header and * C the names of the input files. This program is designed to handle * C MRC format maps but could be modified to handle any map format. * C Map size is currently limited to 16384 points. * C * C Version 1.0 8-AUG-85 KAT for VAX * C Version 1.1 15-AUG-85 TAC + output of average map * C VERSION 2.0 15-APR-86 RAM NEW FORMULA C******************************************************************************** C PARAMETER (MAP$ = 16384) DIMENSION DENS(MAP$),SDENS(MAP$),SDENS2(MAP$),VARIANCE(MAP$) EQUIVALENCE ( DENS , VARIANCE ) INTEGER*4 NAVG,NIN,NVAR,NX,NY,NZ,NXYZ(3),MXYZ(3),MODE CHARACTER IFILE*60,THEDATE*9,THETIME*9 REAL*4 LABELOUT(20,10,2),LABELIN(20,10),TITLE(20) LOGICAL IFIRST EQUIVALENCE (NX,NXYZ(1)), (NY,NXYZ(2)), (NZ,NXYZ(3)) C DATA IFIRST/.TRUE./,NVAR/1/,NIN/2/,NAVG/3/ C CALL TIME(THETIME) !Record date and time of run CALL DATE(THEDATE) WRITE(6,1000) !Print output header READ(5,1002) TITLE !Read title for output map CALL IMOPEN(NVAR,'VAR','NEW') !Open variance map file CALL IMOPEN(NAVG,'AVG','NEW') !Open average file NMAPS = 0 !Zero map counter 1 READ(5,1005,ERR=10,END=100) IFILE !Read input file name GO TO 20 10 WRITE(6,1001) IFILE !Branch to hear for error on input ile name GO TO 1 20 CALL IMOPEN(NIN,IFILE,'RO') !Open input map file CALL IRDHDR(NIN,NXYZ,MXYZ,MODE,DMIN,DMAX,DMEAN) !Read map header IF (IFIRST) THEN !Here we set up header for output maps IF( ( NX * NY ) .GT. MAP$ ) THEN !Will map fit into array dimension WRITE(6,1003) NX,NY,MAP$ !Error message GO TO 200 !EXIT if maps don't fit END IF CALL ITRHDR(NVAR,NIN) !Transfer header input map to output map CALL ITRHDR(NAVG,NIN) !Transfer header to scratch file ENCODE(80,1010,LABELOUT(1,2,1)) THEDATE,THETIME !Create second map title ENCODE(80,1011,LABELOUT(1,2,2)) THEDATE,THETIME !Create second map title DO I=1,20 LABELOUT(I,1,1) = TITLE(I) !transfer title to output map title LABELOUT(I,1,2) = TITLE(I) !transfer title to output map title END DO NNXX = NX !Save dimensions of map to NNYY = NY ! compare with following NNZZ = NZ ! maps NPTS = NX*NY NLABELS = 2 DO ISEC = 1 , NZ !Zero all densities in output CALL IWRSEC(NVAR,SDENS) ! and AVG CALL IWRSEC(NAVG,SDENS2) END DO IFIRST = .FALSE. END IF C IF (NNXX .NE. NX .AND. NNYY .NE. NY .AND. NNZZ .NE. NZ) THEN WRITE(6,1004) NNXX,NNYY,NNZZ !Check to see if map dimensions match GO TO 1 END IF C CALL IRTLAB(NIN,LABELIN,NL) !Transfer first title from each NLABELS = NLABELS + 1 ! input map to the output IF(NLABELS.LE.10) THEN ! variance map. Up to 8 DO JJ=1,20 ! are transferred to the title LABELOUT(JJ,NLABELS,1) = LABELIN(JJ,1) ! field of output variance map LABELOUT(JJ,NLABELS,2) = LABELIN(JJ,1) ! field of output variance map END DO END IF C IY = 0 !Pointer for IY, always zero NMAPS = NMAPS + 1 !Increment map counter DO I = 1 , NNZZ !Loop through sections ISEC = I - 1 !Section pointer always starts at 0 CALL IMPOSN(NIN,ISEC,IY) !Position for reading section ISEC CALL IRDSEC(NIN,DENS,*99) !Read section CALL IMPOSN(NVAR,ISEC,IY) !Position for reading section ISEC CALL IRDSEC(NVAR,SDENS,*99) !Read section CALL IMPOSN(NAVG,ISEC,IY) !Position for reading section ISEC CALL IRDSEC(NAVG,SDENS2,*99) !Read section DO J = 1 , NPTS !Loop through densities SDENS(J) = SDENS(J) + DENS(J) !Summation of densities SDENS2(J) = SDENS2(J) + DENS(J)**2 !Summation of densities squared END DO CALL IMPOSN(NVAR,ISEC,IY) !Position for writing section ISEC CALL IWRSEC(NVAR,SDENS) !Write section CALL IMPOSN(NAVG,ISEC,IY) !Position for writing section ISEC CALL IWRSEC(NAVG,SDENS2) !Write section END DO CALL IMCLOSE(NIN) !Close this input map GO TO 1 100 DMEAN = 0.0 DMAX = -1E7 DMIN = 1E7 IDEG = NMAPS - 1 AVNUM =FLOAT(NMAPS)/FLOAT(IDEG) WRITE(6,1012) NMAPS,IDEG,AVNUM 1012 FORMAT(2I10,F10.4) DO NSEC = 1 , NNZZ !Loop through sections ISEC = NSEC -1 !Section pointer always starts at zero CALL IMPOSN(NVAR,ISEC,IY) !Position for reading section ISEC CALL IRDSEC(NVAR,SDENS,*99) !Read section CALL IMPOSN(NAVG,ISEC,IY) !Position for reading section ISEC CALL IRDSEC(NAVG,SDENS2,*99) !Read section DO J = 1 , NPTS !Calculate variance point by point TMP1=SDENS2(J)/IDEG TMP2=(SDENS(J)/NMAPS ) **2 VARIANCE(J) = TMP1 - (AVNUM * TMP2) IF ( VARIANCE(J) .GT. 0.0 ) THEN VARIANCE(J) = SQRT(VARIANCE(J)) ELSE c if (variance(j) .lt. 0.0) write(6,1013) VARIANCE(J) = 0.0 END IF IF ( VARIANCE(J) .LT. DMIN ) DMIN = VARIANCE(J)!Find limits IF ( VARIANCE(J) .GT. DMAX ) DMAX = VARIANCE(J)!Find limits DMEAN = DMEAN + VARIANCE(J) !Find mean C SDENS(J)=SDENS(J)/NMAPS !AVERAGE MAP IF ( SDENS(J) .LT. AMIN ) AMIN = SDENS(J) !Find limits IF ( SDENS(J) .GT. AMAX ) AMAX = SDENS(J) !Find limits AMEAN = AMEAN + SDENS(J) !Find mean C END DO C WRITE(6,1012) NMAPS,IDEG,AVNUM C1012 FORMAT(2I10,F10.4) CALL IMPOSN(NVAR,ISEC,IY) !Position for writing section CALL IWRSEC(NVAR,VARIANCE) !Write section C CALL IMPOSN(NAVG,ISEC,IY) !Position for writing section CALL IWRSEC(NAVG,SDENS) !Write section END DO C CLOSE (UNIT=NSCRATCH,STATUS='DELETE') !Delete AVG file DMEAN = DMEAN / NPTS / NNZZ !Calculate mean AMEAN = AMEAN / NPTS / NNZZ !Calculate mean C CALL IALLAB(NVAR,LABELOUT(1,1,1),NLABELS) !Transfer labels to output map CALL IWRHDR(NVAR,TITLE,-1,DMIN,DMAX,DMEAN) !Write output map header CALL IMCLOSE(NVAR) !Close output variance map C CALL IALLAB(NAVG,LABELOUT(1,1,2),NLABELS) !Transfer labels to output map CALL IWRHDR(NAVG,TITLE,-1,AMIN,AMAX,AMEAN) !Write output map header CALL IMCLOSE(NAVG) !Close output AVG GO TO 200 C 99 WRITE(6,1006) C 1000 FORMAT(/,' VARIMAP, Version 1.1: Calculate density variance ', $'and average between two or more maps') 1001 FORMAT(/,' !!! CANNOT OPEN INPUT FILE: ',A/) 1002 FORMAT(20A4) 1003 FORMAT(//,' Map dimensions too large for program, NX =',I5, A ', NY =',I5,', map limit =',I10) 1004 FORMAT(' !! Dimensions of input map do not match output map!!') 1005 FORMAT(A) 1006 FORMAT(/' END OF FILE ON MAP INPUT FILE !!!!') 1010 FORMAT('VARIMAP, Version 1.1: Calculate variance map',A,', ',A) 1011 FORMAT('VARIMAP, Version 1.1: Calculate average map ',A,', ',A) 1013 format(' negative variance encountered !!!!!') C 200 CALL EXIT END