Remember: This thesis is from the early eighties. Computers were quite different then.
Gridding is the interpolation of irregularly spaced data into an evenly spaced matrix. Many methods have been developed to perform gridding, including various types of weighted averages (such as Davis, 1973), linear fitting (Palmer, 1969, and Walters, 1969), cubic splines (Bhattacharyya, 1969), piecewise continuous cubic polynomials (Akima, 1974, Hessing and others, 1972), polynomial surface fitting (Braile, 1978), and difference equations for minimum curvature (Swain, 1976). Whether or not a gridding method is good or not depends on the following factors: how well short wavelength features are represented; the accuracy in regions of sparse data; the interpolation of the data near the edges of the grid; and how much computer time is required (Braile, 1978). Each of these methods can give good results, but it often depends on the amount and distribution of the irregularly spaced data. Braile (1978) compares many of these gridding methods and determines that polynomial surface fitting is the best method of gridding.
The program presented here uses the polynomial surface fitting method in Braile (1978). Polynomial surface fitting uses data points within a region to calculate the coefficients of an Nth order polynomial to determine the value of grid points located in the center of the region. The number of data points within the region required to calculate the polynomial depends on the desired order of the polynomial. The coefficients to a single global polynomial, one that covers the entire data set, could be calculated so that the value at each grid location could be calculated. However, this would tend to filter out much of the fine detail in the data set. It is more useful to caculate many regional polynomials, with each covering only a small portion of the data set. To calculate the values of a few grid locations, the data set would be searched for all data points within a radius R of the grid locations. Only these would be used to calculate the coefficients of the polynomial. The values at the grid locations would be calculated from these coefficients. This process would be repeated for each of the grid locations. A significant problem with this is the amount of time spent searching the data set for the nearest data points. For every regional polynomial, the entire data set would need to be searched for the closest data points. However, a more significant problem is the distribution of the data points within the polynomial region. If within the radius R the data points are not well distributed, the polynomial can become unstable and shoot to + or - infinity.
The program presented here avoids both of these problems very well. To create a matrix of N row by M columns, the data is sorted into a very large 'pegboard' of 5N rows by 5M columns. For each data point's X, Y location, the data is placed into the nearest position in the pegboard. Each position in the pegboard contains 3 variables to store the actual X, Y location and the Z gravity value. If more than one data value is to be placed in a position in the pegboard, the values are averaged together. This insures the data is not concentrated within a small area.
Once the data is sorted into the pegboard, the grid can be calculated. For each grid location, all data within a radius R are needed to calculate the polynomial coefficients. All pegboard positions withing the radius R are checked for data. Many of the positions will be empty, but since only a limited region of the data set needs to be searched, it is very quick. If not enough data exists within this radius R, R can be increased. Once the coefficients of the polynomial have been determined, the value at the grid location can be calculated. If the value does not seem correct, (such as not within a minimum and maximum acceptable range), the program can alter the radius or the order of the polynomial and try again.
The program GRID_USING_POLY_FIT requires approximately 2.5 Mb of storage. Obviously it cannot hold all of its arrays in memory at one time. It is because of the virtual memory characteristics of the VAX and similar computers that make this method so efficient. With virtual memory, the arrays are stored on a work space on the system disk, and the computer reads and writes the pages of memory as needed. This is completely transparrent to the user.
FORTRAN-77 Computer Program GRID_USING_POLY_FIT
PROGRAM GRID_USING_POLY_FIT
REAL MAP(3,560,320),GRID(128),XYZ(3),XY(3,400),XC(137)
CHARACTER FNAME*32,FNAME1*32,QUES
INTEGER*2 NCOUNT(560,320)
INTEGER*4 NCOEF(15),IRANGE(4)
LOGICAL*1 ERROR,GOOD_ANS
DATA NCOEF/3,6,10,15,21,28,36,45,55,66,78,91,105,120,136/
DATA IRANGE/10,20,30,40/
C
TYPE *,' Name of gridded XYZ map file?'
ACCEPT 10,FNAME
10 FORMAT(A32)
TYPE *,' Name for NEW map file?'
ACCEPT 10,FNAME1
TYPE *,' Coordinate of lower left corner of input data?'
ACCEPT *,XMIN,YMIN
TYPE *,' Coordinate of lower left corner of output matrix?'
ACCEPT *,X0,Y0
TYPE *,' INPUT DATA: Number of cells per kilometer?'
ACCEPT *,DELTA
TYPE *,' OUTPUT MATRIX: Number of columns?'
ACCEPT *,NX
TYPE *,' OUTPUT MATRIX: Number of rows?'
ACCEPT *,NY
TYPE *,' OUTPUT MATRIX: Distance between rows & columns?'
ACCEPT *,DELXY
TYPE *,' Minimum and maximum acceptable answers?'
ACCEPT *,ANSMIN,ANSMAX
C
C INITIALIZE MATRICES
C
DO J=1,320
DO I=1,560
NCOUNT(I,J)=0
DO L=1,3
MAP(L,I,J)=0
END DO
END DO
END DO
C
C Load data matrix
C
OPEN(UNIT=1,FILE=FNAME,STATUS='OLD',
> ACCESS='DIRECT',RECORDSIZE=3)
READ(1'1) NUM
IREC=1
ERROR=.TRUE.
DO IR=1,NUM-1
IREC=IREC+1
READ(1'IREC,ERR=89) (XYZ(N),N=1,3)
I=(XYZ(1)-XMIN)*DELTA+1.
J=(XYZ(2)-YMIN)*DELTA+1.
IF(I.GT.0.AND.J.GT.0.AND.I.LE.560.AND.J.LE.320) THEN
NCOUNT(I,J)=NCOUNT(I,J)+1
DO N=1,3
MAP(N,I,J)=MAP(N,I,J)+XYZ(N)
END DO
ENDIF
END DO
ERROR=.FALSE.
89 IF(ERROR) TYPE *,' READ ERROR AT RECORD',IREC
CLOSE(UNIT=1)
DO J=1,320
DO I=1,560
IF(NCOUNT(I,J).GT.1) THEN
RN=NCOUNT(I,J)
DO N=1,3
MAP(N,I,J)=MAP(N,I,J)/RN
END DO
ENDIF
END DO
END DO
C
C Carry out polynomial fitting
C
OPEN(UNIT=1,FILE=FNAME1,STATUS='NEW',
> ACCESS='DIRECT',RECORDSIZE=NX)
TYPE *,' BEGIN POLYNOMIAL FITTING'
DO IY=1,NY
Y1=Y0+FLOAT(IY-1)*DELXY
IY1=(Y1-YMIN)*DELTA+1.
DO IX=1,NX
X1=X0+FLOAT(IX-1)*DELXY
IX1=(X1-XMIN)*DELTA+1.
ITRY=0
GOOD_ANS=.FALSE.
DO WHILE(.NOT.GOOD_ANS.AND.ITRY.LT.3)
ITRY=ITRY+1
NUM=0
DO LY=-IRANGE(ITRY)+1,IRANGE(ITRY)
J=IY1+LY
IF(J.GT.0.AND.J.LE.320) THEN
DO LX=-IRANGE(ITRY)+1,IRANGE(ITRY)
I=IX1+LX
IF(I.GT.0.AND.I.LE.560.AND.NUM.LT.400) THEN
IF(NCOUNT(I,J).GT.0) THEN
NUM=NUM+1
XY(1,NUM)=MAP(1,I,J)-X0
XY(2,NUM)=MAP(2,I,J)-Y0
XY(3,NUM)=MAP(3,I,J)
ENDIF
ENDIF
END DO
ENDIF
END DO
C
ERROR=.TRUE.
NORD=5
RNUM=NUM
DO WHILE(FLOAT(NCOEF(NORD))*1.3.GT.RNUM)
NORD=NORD-1
IF(NORD.LT.2) GOTO 999
END DO
ERROR=.FALSE.
CALL LSPOLY(XY,XC,NUM,NORD)
ANS=POLY(XC,X1-X0,Y1-Y0,NORD)
999 IF(ERROR) TYPE *,'ERROR -- Not enought points ',IX,IY
GOOD_ANS=ANS.GT.ANSMIN.AND.ANS.LT.ANSMAX
END DO
IF(.NOT.GOOD_ANS) TYPE *,' Bad answer after 4 tries'
GRID(IX)=ANS
END DO
C
C WRITE OUTPUT FILE
C
WRITE(1'IY) (GRID(IX),IX=1,NX)
END DO
CLOSE(UNIT=1)
STOP
END
SUBROUTINE LSPOLY(XY,XC,NSTA,NORD1)
REAL XY(3,400),A(400,137),XC(137),B(400)
C****************************************************************
C LSQSRF
C****************************************************************
C
C LSQSRF computes a least squares polynomial surface of
C order NORD1 to data given as B(X,Y) where NSTA data
C points are input. Terms in polynomial are as follows:
C
C NORD NCOEF TERMS
C ---- ----- -----
C 1 3 Z = const + a X + b Y
C 2 6 + c X**2 + d X*Y + e Y**2
C 3 10 + f X**3 + g X**2*Y+ h X*Y**2 + i Y**3
C 4 15 + j X**4 + k X**3*Y+ l X**2*Y**2
C + m X*Y**3 + n Y**4
C ETC.
C
NCOEF=FLOAT(((NORD1+1)*(NORD1+2)/2))
DO L1=1,NSTA
NV=2
X=XY(1,L1)
Y=XY(2,L1)
B(L1)=XY(3,L1)
A(L1,1)=X
A(L1,2)=Y
A(L1,NCOEF)=1.
DO M1=2,NORD1
DO K=1,M1
NV=NV+1
L=NV-M1
A(L1,NV)=A(L1,L)*X
END DO
NV=NV+1
A(L1,NV)=Y**M1
END DO
END DO
C
C Solve for A(m,n)*XC(n)=B(m)
C
CALL GOLUB(A,XC,B,NSTA,NCOEF)
RETURN
END
SUBROUTINE GOLUB(A,X,B,M,N)
C
C A(M,N) ; B(M) given with M>N solves for X(N) such that
C || B - AX || = MINIMUM
C
C Method of G.GOLUB, NUMERSICHE MATHEMATIK 7, 206-216 1965
C from CLAREBOUT p. 119
C
C dimension U(M)
C perform N orthogonal transformations to A(.,.) to
C upper triangularize the matrix
C
REAL A(400,137),X(137),B(400),U(400)
DOUBLE PRECISION DSUM,DAJ,DAI,DSIGMA,DBI,D1
D1=1.0
DO K=1,N
DSUM=0.0
DO I=K,M
DAJ=A(I,K)
DSUM=DSUM+DAJ*DAJ
END DO
DAI=A(K,K)
DSIGMA=DSIGN(DSQRT(DSUM),DAI)
DBI=DSQRT(D1+DAI/DSIGMA)
FACT=D1/(DSIGMA*DBI)
U(K)=DBI
DO I=K+1,M
U(I)=FACT*A(I,K)
END DO
C
C I - UU' is a symmetric, orthogonal matrix which when applied
C to A(.,.) will zero the elements below the diagonal K
C
C Apply the orthoganal transformation
C
DO J=K,N
FACT=0.0
DO I=K,M
FACT=FACT+U(I)*A(I,J)
END DO
DO I=K,M
A(I,J)=A(I,J)-FACT*U(I)
END DO
END DO
FACT=0.0
DO I=K,M
FACT=FACT+U(I)*B(I)
END DO
DO I=K,M
B(I)=B(I)-FACT*U(I)
END DO
END DO
C
C Back substitute to recursively yield X(.)
C
X(N)=B(N)/A(N,N)
LIM=N-1
DO I=1,LIM
IROW=N-I
SUM=0.0
DO J=1,I
SUM=SUM+X(N-J+1)*A(IROW,N-J+1)
END DO
X(IROW)=(B(IROW)-SUM)/A(IROW,IROW)
END DO
RETURN
END
REAL FUNCTION POLY(XC,X,Y,NORD1)
C calculate value of least squares
C polynomial surface at X,Y
REAL C(137),XC(137)
NCOEF=FLOAT(((NORD1+1)*(NORD1+2)/2))
C(1)=X
C(2)=Y
C(NCOEF)=1.
NV=2
DO M1=2,NORD1
DO K=1,M1
NV=NV+1
L=NV-M1
C(NV)=C(L)*X
END DO
NV=NV+1
C(NV)=Y**M1
END DO
SUM=0.0
DO K=1,NCOEF
SUM=SUM+C(K)*XC(K)
END DO
POLY=SUM
RETURN
END
54 69 6D 20 46 6F 67 61 72 74 79 20 6C 6F 76 65 73 20
53 74 65 76 65 20 52 69 6F 72 64 61 6E 21
[Editors note: As it is bounded and sitting on some shelf in the Geoscience Library, I placed the above coded message in various locations in my thesis. As I was 25 and only recently out, this was making quite a statement. $10 or a big kiss, which every you prefer, to the first person who emails me the correct uncoded message.]
Since different types of wavenumber filtering is usually performed on a single data set, many of the steps would be performed over and over. To eliminate duplication of some of the calculations, all steps up to creation of the ideal filter matrix is placed in the program SETUP. This program would be used once and the results given to the program FILTER.
This program, written in VAX/VMS Fortran-77, performs the following steps:
The arrays of the program are currently dimensioned to 256 by 128. Since the data set is to be expanded to limit edge effects, the maximum size of the input data matrix is 128 columns by 64 rows. Note that the Fast Fourier Transform algorithm listed here requires the number of rows and columns be powers of 2.
For efficient I/O, the input data matrix must be stored as a direct access file of recordsize equal to the number of columns. The following VAX/VMS Fortran-77 code will produce such a file.
REAL MATRIX(128,64)
CHARACTER FNAME*32
.
.
.
TYPE *,' Name of output file?'
ACCEPT 10,FNAME
10 FORMAT(A32)
TYPE *,' Number of columns?'
ACCEPT *,NCOL
TYPE *,' Number of rows?'
ACCEPT *,NROW
.
.
.
OPEN(UNIT=1,FILE=FNAME,STATUS='NEW',ACCESS='DIRECT',
> RECORDSIZE=NCOL)
DO NR=1,NROW
WRITE(1'NR) (MATRIX(NC,NR),NC=1,NCOL)
END DO
CLOSE(UNIT=1)
.
.
.
In other versions of Fortran-77, one must specify the recordsize as number of bytes rather than number of words. Such code would be RECORDSIZE=NCOL*4 rather than RECORDSIZE=NCOL. The WRITE statement could also be of the form WRITE(1,REC=NR)...
The program will ask the user to input the following.
PROGRAM SETUP_FOR_FILTER
C
C This program sets up the grid file for filter.
C It calculates the FFT of the spatial data.
C
CHARACTER*32 IN,OUT
COMMON/BLK/D1(2,256,128),NX,NY,DELTA
T1=SECNDS(0.0) ! VAX RTL ROUTINE
C
C Read in parameters.
C
TYPE *,' Set up matrix file for filtering -- DAF = NX'
TYPE *,' NAME OF DIRECT ACCESS INPUT FILE?'
ACCEPT 10,IN
10 FORMAT(A32)
TYPE *,' NAME OF DIRECT ACCESS OUTPUT FILE?'
ACCEPT 10,OUT
TYPE *,' NUMBER OF COLUMNS?'
ACCEPT *,NX2
TYPE *,' NUMBER OF ROWS?'
ACCEPT *,NY2
TYPE *,' GRID INTERVAL?'
ACCEPT *,DELTA
TYPE *,' DO YOU WISH TO EXTEND MATRIX BY'
TYPE *,' 1...EXTENDING END VALUES'
TYPE *,' 2...REFLECTING TO FOUR QUADRANTS'
TYPE *,'?'
ACCEPT *,IR
NX=NX2+NX2
RNX=NX
NY=NY2+NY2
RNY=NY
OPEN(UNIT=1,NAME=IN,TYPE='OLD',ACCESS='DIRECT',
> READONLY,RECORDSIZE=NX2)
DO IY=1,NY2
READ(1'IY) (D1(1,IX,IY),IX=1,NX2)
END DO
CLOSE(UNIT=1)
C----------------------------------------------------------------
C REFLECT OR EXTEND LAST VALUES
C----------------------------------------------------------------
IF(IR.EQ.2) THEN
CALL REFLECT
ELSE
CALL EXTEND
ENDIF
C----------------------------------------------------------------
C CALCULATE AND REMOVE MEAN VALUE OF DATA MATRIX
C----------------------------------------------------------------
CALL REMOVE(XMEAN)
TYPE *,' MEAN WAS ',XMEAN
C----------------------------------------------------------------
C COMPUTE TRANSFORM OF DATA
C----------------------------------------------------------------
CALL FFT2D(-1)
C----------------------------------------------------------------
C CREATE OUTPUT FILE
C----------------------------------------------------------------
IREC_LEN=NX2*4
OPEN(UNIT=2,NAME=OUT,TYPE='NEW',ACCESS='DIRECT',
> RECORDSIZE=IREC_LEN)
WRITE(2'1) XMEAN,DELTA,RNX,RNY,D1(1,1,1)
DO J=1,NY
WRITE(2'J+1) (D1(1,I,J),D1(2,I,J),I=1,NX)
END DO
CLOSE(UNIT=2)
TYPE *,' FFT ERROR WAS',D1(1,1,1)
ELAPSE=SECNDS(T1) ! VAX RTL ROUTINE
TYPE *,' ELAPSE TIME',ELAPSE
STOP
END
SUBROUTINE FFT2D(NSIGN)
COMPLEX C(256),H,CTEMP
COMMON/BLK/H(256,128),NX,NY,DELTA
C
C INCORE METHOD
C FFT2D PERFORMS BOTH A FORWARD OR INVERS FAST FOURIER
C TRANSFORM. FFT2D IS THE DRIVER THAT PASSES THE CORRECT
C VECTORS TO FORK WITH PERFORMS THE ACTUAL TRANSFORMING
C
C NX NUMBER OF COLUMNS
C NY NUMBER OF ROWS
C NSIGN DIRECTION OF DESIRED TRANSFORMATION
C +1 FREQUENCY TO SPATIAL
C -1 SPACIAL TO FREQUENCY
C
SIGNI=FLOAT(NSIGN)
C
C OPERATE BY ROWS
C
DO IY=1,NY
DO IX=1,NX
C(IX)=H(IX,IY)
END DO
CALL FORK(C,NX,SIGNI)
DO IX=1,NX
H(IX,IY)=C(IX)
END DO
END DO
IF(NY.EQ.1) RETURN
C
C OPERATE BY COLUMNS
C
DO IX=1,NX
DO IY=1,NY
C(IY)=H(IX,IY)
END DO
CALL FORK(C,NY,SIGNI)
DO IY=1,NY
H(IX,IY)=C(IY)
END DO
END DO
RETURN
END
SUBROUTINE FORK(CX,LXX,SIGNI)
COMPLEX CX(256),CW,CTEMP,CON2
C
C FAST FOURIER TRANSFORM, MODIFIED FROM CLAERBOUT, J.F.,
C FUNDAMENTALS OF GEOPHYSICAL DATA PROCESSING,
C MCGRAW-HILL 1976
C
C CX DATA VECTOR TO BE TRANSFORMED
C LXX LENGTH OF DATA VECTOR CX TO BE TRANSFORMED
C MUST BE A POWER OF 2 (LXX=2**INTEGER)
C SIGNI DIRECTION OF TRANSFORMATION
C +1 INVERSE--FREQUENCY TO SPATIAL
C -1 FORWARD--SPATIAL TO FREQUENCY
C
C NORMALIZATION PERFORMED BY DIVIDING BY DATA LENGTH
C UPON THE FORWARD TRANS.
C
LX=LXX
LXH=LX/2
J=1
DO 10 I=1,LX
IF(I.LT.J) THEN
CTEMP=CX(J)
CX(J)=CX(I)
CX(I)=CTEMP
ENDIF
M=LXH
DO WHILE(M.GE.1)
IF(J.LE.M) GOTO 10
J=J-M
M=M/2
END DO
10 J=J+M
L=1
DO WHILE(L.LT.LX)
ISTEP=2*L
CON2=(0.0,3.14159265)/FLOAT(L)*SIGNI
DO M=1,L
CW=CEXP(CON2*FLOAT(M-1))
DO I=M,LX,ISTEP
CTEMP=CW*CX(I+L)
CX(I+L)=CX(I)-CTEMP
CX(I)=CX(I)+CTEMP
END DO
END DO
L=ISTEP
END DO
IF(SIGNI.LE.0.0) THEN
SC=1./FLOAT(LX)
DO I=1,LX
CX(I)=CX(I)*SC
END DO
ENDIF
RETURN
END
SUBROUTINE EXTEND
COMMON/BLK/D1(2,256,128),NXD,NYD,DELTA
C
C THIS ROUTINE EXTENDS THE END VALUES
C
C
C ____________________
C | | | |
C |2 |3 |1 |
C ____________________ --------------------
C | | | | | | | |
C |7 |6 |5 | |6 |5 |7 |
C -------------------- --------------------
C | |or | | axis |or | | |
C | | ig | | shifted | ig | | |
C | | in | | from | in | | |
C |8 | al|4 | <=== | al|4 |8 |
C -------------------- to --------------------
C | | | | ===>
C |1 |2 |3 |
C --------------------
C
C
C
C
NX=NXD/2
NY=NYD/2
NX2=NX/2
NY2=NY/2
DO I=1,NX
R1=D1(1,I,1)
R2=D1(2,I,NY)
DO J=1,NY2
D1(1,I,NY+J)=R2
D1(1,I,NY+NY2+J)=R1
END DO
END DO
DO J=1,NX
R1=D1(1,1,J)
R2=D1(1,NX,J)
DO I=1,NX2
D1(1,NX+I,J)=R2
D1(1,NX+I+NX2,J)=R1
END DO
END DO
R1=D1(1,1,1)
R2=D1(1,NX,1)
R3=D1(1,NX,NY)
R4=D1(1,1,NY)
DO I=1,NX2
IX=NX+I
IX2=IX+NX2
DO J=1,NY2
JY=NY+J
JY2=JY+NY2
D1(1,IX,JY)=R3
D1(1,IX,JY2)=R2
D1(1,IX2,JY)=R4
D1(1,IX2,JY2)=R1
END DO
END DO
RETURN
END
SUBROUTINE REFLECT
COMMON/BLK/D1(2,256,128),NX,NY,DELTA
C
C THIS ROUTINE REFLECTS THE DATA IN ONE QUADRANT
C INTO THE OTHER THREE
C
NX2=NX/2
NY2=NY/2
DO I=1,NX2
IX=NX-I+1
DO J=1,NY2
JY=NY-J+1
TEMP=D1(1,I,J)
D1(1,IX,J)=TEMP
D1(1,IX,JY)=TEMP
D1(1,I,JY)=TEMP
END DO
END DO
RETURN
END
C
SUBROUTINE REMOVE(XMEAN)
REAL*8 SUM,RNUM,DD
COMMON/BLK/D1(2,256,128),NX,NY,DELTA
SUM=0.0
DO I=1,NX
DO J=1,NY
DD=D1(1,I,J)
SUM=SUM+DD
END DO
END DO
RNUM=NX*NY
XMEAN=SUM/RNUM
DO I=1,NX
DO J=1,NY
D1(1,I,J)=D1(1,I,J)-XMEAN
D1(2,I,J)=0.0
END DO
END DO
RETURN
END
54 69 6D 20 46 6F 67 61 72 74 79 20 6C 6F 76 65 73 20
53 74 65 76 65 20 52 69 6F 72 64 61 6E 21
[Editors note: As it is bounded and sitting on some shelf in the Geoscience Library, I placed the above coded message in various locations in my thesis. As I was 25 and only recently out, this was making quite a statement. $10 or a big kiss, which every you prefer, to the first person who emails me the correct uncoded message.]
This program performs various types of filtering on a matrix that has been prepared by SETUP. This program, also written in VAX/VMS Fortran-77, performs the following steps:
Four types of filters are available. They are Band-pass, Strike-pass, Upward/Downward continuation, and Nth order vertical or horizontal derivative.
The program will ask the user to input the following.
If band-pass filtering is selected, the program will ask for:
If strike-pass filtering is selected, the program will ask:
If Nth order derivative filtering is selected:
If upward or downward continuation filtering is selected:
For all filters, the program will be asked for:
This program filters a data set of 128 by 64 extended to 256 by 128 by SETUP in less than one minute of CPU time on the VAX 11-750.
If the user wishes to perform multiple filters on a data set, such as a first order XY derivative or a downward continuation of a band-passed data set, the user must perform the first filter, pass the results through the program SETUP, and then perform the second filter.
PROGRAM FILTER
C
C THIS PROGRAM RUNS THE FILTERING ROUTINES IN THE WAVENUMBER
C DOMAIN OF UNIFORMLY GRIDDED ONE- OR TWO-DIMENSIONAL DATA SETS
C
C FILTERS ARE BAND-PASS
C NTH DERIVATIVE
C STRIKE PASS/REJECT
C UPWARD/DOWNWARD CONTINUATION
C
DIMENSION D1(2,256,128)
COMPLEX H,CSTORE(256)
CHARACTER DBUF*9,TBUF*9,QUES
CHARACTER*32 IN,OUT,LST
LOGICAL*2 ADM
EQUIVALENCE (D1,H)
COMMON/BLK/H(256,128),NX,NY,DELTA
T1=SECNDS(0.0)
CALL DATE(DBUF)
CALL TIME(TBUF)
C
C READ IN PARAMETERS.
C
TYPE *,' Name of file from program SETUP?'
ACCEPT 10,IN
10 FORMAT(A32)
TYPE *,' Name of output file?'
ACCEPT 10,OUT
IS=INDEX(OUT,'.') ! ASSUMES NO SUBDIRECTORY INCLUDED
LST=OUT(1:IS)//'FIL'
TYPE 30,LST
30 FORMAT(' List file is ',A32)
OPEN(UNIT=4,NAME=LST,STATUS='NEW')
WRITE(4,40) DBUF,TBUF
40 FORMAT('0 IN CORE FILTERING PROGRAM',10X,A9,10X,A9)
WRITE(4,50) IN,OUT
50 FORMAT(' INPUT FILE: ',A32,10X,'OUTPUT FILE: ',A32)
TYPE *,' Number of columns in original matrix ?'
ACCEPT *,NX2
IREC_LEN=NX2*4
C
TYPE *,' Which filter do you wish to perform?'
TYPE *,' '
TYPE *,' 1 BAND-PASS'
TYPE *,' 2 STRIKE-PASS'
TYPE *,' 3 Nth DERIVATIVE'
TYPE *,' 4 UPWARD/DOWNWARD CONTINUATION'
TYPE *,'?'
ACCEPT *,IFIL
TYPE *,' '
OPEN(UNIT=2,NAME=IN,STATUS='OLD',ACCESS='DIRECT',
> RECORDSIZE=IREC_LEN,READONLY)
READ(2'1) XMEAN,DELTA,RNX,RNY,FFTERR
C****************************************************************
IF(IFIL.EQ.1) THEN
TYPE *,' BAND-PASS/REJECT'
TYPE *,' '
TYPE *,' Range of wavenumbers is 0.0 to',.5/DELTA
TYPE *,' '
TYPE *,' Range of wavenumbers to be passed or rejected?'
ACCEPT *,CUTLO,CUTHI
IF(CUTHI.LT.CUTLO) THEN
TEMP=CUTHI
CUTHI=CUTLO
CUTLO=TEMP
ENDIF
TYPE *,' DO YOU WISH TO'
TYPE *,' -1 Reject wavenumbers between hi and low'
TYPE *,' 1 Pass wavenumbers between hi and low'
TYPE *,'?'
ACCEPT *,NPASS
WRITE(4,*) ' BAND-PASS FILTER'
IF(NPASS.EQ.1) THEN
WRITE(4,*) ' PASS WAVENUMBERS BETWEEN',CUTLO,CUTHI
ELSE
WRITE(4,*) ' REJECT WAVENUMBERS BETWEEN',CUTLO,CUTHI
ENDIF
C****************************************************************
ELSEIF(IFIL.EQ.2) THEN
TYPE *,' STRIKE PASS/REJECT'
TYPE *,' '
TYPE *,' Smallest angle?'
ACCEPT *,ANG1
TYPE *,' Largest angle?'
ACCEPT *,ANG2
TYPE *,' DO YOU WISH TO'
TYPE *,' -1 Reject azimuths between angles'
TYPE *,' 1 Pass azimuths between angles'
TYPE *,'?'
ACCEPT *,NPASS
WRITE(4,*) ' STRIKE-PASS FILTER'
IF(NPASS.EQ.1) THEN
WRITE(4,*) ' PASS ANGLES BETWEEN',ANG1,ANG2
ELSE
WRITE(4,*) ' REJECT ANGLES BETWEEN',ANG1,ANG2
ENDIF
C****************************************************************
ELSEIF(IFIL.EQ.3) THEN
TYPE *,' Nth DERIVATIVE'
TYPE *,' '
TYPE *,' Order of spatial derivative?'
ACCEPT *,NTH
TYPE *,' DIRECTION IN WHICH TO CALCULATE THE DERIVATIVE'
TYPE *,' 1 Vertical derivative'
TYPE *,' 2 Horizontal derivative in X direction'
TYPE *,' 3 Horizontal deirvative in Y direction'
TYPE *,'?'
ACCEPT *,NWAY
IF(NWAY.EQ.1) THEN
WRITE(4,41) NTH,'VERTICAL'
ELSEIF(NWAY.EQ.2) THEN
WRITE(4,41) NTH,'HORIZ. X'
ELSE
WRITE(4,41) NTH,'HORIZ. Y'
ENDIF
41 FORMAT(I3,' ORDER DERIVATIVE FILTER IN ',A8,' DIRECTION')
C****************************************************************
ELSE
TYPE *,' UPWARD/DOWNWARD CONTINUATION'
TYPE *,' '
TYPE *,' Distance UPWARD to continue the data?'
ACCEPT *,ZCON
WRITE(4,*) ' UP/DOWN CONTINUED TO',ZCON
ZCON=-ZCON
C****************************************************************
ENDIF
TYPE *,' Do you wish the mean of the data to be added'
TYPE *,' back for the final results?'
ACCEPT 60,QUES
60 FORMAT(1A1)
ADM=QUES.EQ.'Y'.OR.QUES.EQ.'y'
C****************************************************************
C OPEN FILE THAT HAS BEEN PREVIOUSLY SET UP
C****************************************************************
IF(.NOT.ADM) XMEAN=0.0
NX=RNX
NY=RNY
NX2=NX/2
NY2=NY/2
WRITE(4,*) ' EXTENDED ARRAY OF',NX,' COLUMNS BY',NY,' ROWS'
WRITE(4,*) ' ORIGINAL ARRAY OF',NX2,NY2
WRITE(4,*) ' GRID SPACING OF',DELTA,' MEAN OF ',XMEAN
WRITE(4,*) ' FFT ERROR OF ',FFTERR
C----------------------------------------------------------------
C CREATE IDEAL FILTER
C----------------------------------------------------------------
WRITE(4,*)
IF(IFIL.EQ.1) THEN
CALL BNDPAS(CUTLO,CUTHI,NPASS)
ELSEIF(IFIL.EQ.2) THEN
CALL STRIKE(ANG1,ANG2,NPASS)
ELSEIF(IFIL.EQ.3) THEN
CALL DERIV(NTH,NWAY)
ELSE
CALL CONTIN(ZCON)
ENDIF
C----------------------------------------------------------------
C MULTIPLY TRANSFORM OF DATA AND FILTER (CONVOLVING)
C----------------------------------------------------------------
C READ IN EACH VALUE, MULTIPLY
C BY FILTER, STORE RESULTS
IREC=1
DO IY=1,NY
IREC=IREC+1
READ(2'IREC) (CSTORE(IX),IX=1,NX)
DO IX=1,NX
H(IX,IY)=H(IX,IY)*CSTORE(IX)
END DO
END DO
CLOSE(UNIT=2)
C----------------------------------------------------------------
C TRANSFORM RESULTS TO SPATIAL DOMAIN
C----------------------------------------------------------------
CALL FFT2D(1)
C----------------------------------------------------------------
C WRITE OUTPUT FILE (ONLY FIRST QUADRANT)
C----------------------------------------------------------------
OPEN(UNIT=1,NAME=OUT,STATUS='NEW',
> ACCESS='DIRECT',RECORDSIZE=NX2)
DO IY=1,NY2
WRITE(1'IY) (D1(1,IX,IY)+XMEAN,IX=1,NX2)
END DO
C----------------------------------------------------------------
C CLOSE FILES
C----------------------------------------------------------------
CLOSE(UNIT=1)
ELAPSE=SECNDS(T1) ! VAX RTL ROUTINE
WRITE(4,90) ELAPSE
90 FORMAT(' ELAPSE TIME',F14.3)
CLOSE(UNIT=4)
STOP
END
SUBROUTINE FFT2D(NSIGN)
COMPLEX C(256),H,CTEMP
COMMON/BLK/H(256,128),NX,NY,DELTA
C
C INCORE METHOD
C FFT2D PERFORMS BOTH A FORWARD OR INVERS FAST FOURIER
C TRANSFORM. FFT2D IS THE DRIVER THAT PASSES THE CORRECT
C VECTORS TO FORK WITH PERFORMS THE ACTUAL TRANSFORMING
C
C NX NUMBER OF COLUMNS
C NY NUMBER OF ROWS
C NSIGN DIRECTION OF DESIRED TRANSFORMATION
C +1 FREQUENCY TO SPATIAL
C -1 SPACIAL TO FREQUENCY
C
SIGNI=FLOAT(NSIGN)
C
C OPERATE BY ROWS
C
DO IY=1,NY
DO IX=1,NX
C(IX)=H(IX,IY)
END DO
CALL FORK(C,NX,SIGNI)
DO IX=1,NX
H(IX,IY)=C(IX)
END DO
END DO
IF(NY.EQ.1) RETURN
C
C OPERATE BY COLUMNS
C
DO IX=1,NX
DO IY=1,NY
C(IY)=H(IX,IY)
END DO
CALL FORK(C,NY,SIGNI)
DO IY=1,NY
H(IX,IY)=C(IY)
END DO
END DO
RETURN
END
SUBROUTINE FORK(CX,LXX,SIGNI)
COMPLEX CX(256),CW,CTEMP,CON2
C
C FAST FOURIER TRANSFORM, MODIFIED FROM CLAERBOUT, J.F.,
C FUNDAMENTALS OF GEOPHYSICAL DATA PROCESSING,
C MCGRAW-HILL 1976
C
C CX DATA VECTOR TO BE TRANSFORMED
C LXX LENGTH OF DATA VECTOR CX TO BE TRANSFORMED
C MUST BE A POWER OF 2 (LXX=2**INTEGER)
C SIGNI DIRECTION OF TRANSFORMATION
C +1 INVERSE--FREQUENCY TO SPATIAL
C -1 FORWARD--SPATIAL TO FREQUENCY
C
C NORMALIZATION PERFORMED BY DIVIDING BY DATA LENGTH
C UPON THE FORWARD TRANS.
C
LX=LXX
LXH=LX/2
J=1
DO I=1,LX
IF(I.LT.J) THEN
CTEMP=CX(J)
CX(J)=CX(I)
CX(I)=CTEMP
ENDIF
M=LXH
DO WHILE(M.GE.1)
IF(J.LE.M) GOTO 10
J=J-M
M=M/2
END DO
10 J=J+M
END DO
L=1
DO WHILE(L.LT.LX)
ISTEP=2*L
CON2=(0.0,3.14159265)/FLOAT(L)*SIGNI
DO M=1,L
CW=CEXP(CON2*FLOAT(M-1))
DO I=M,LX,ISTEP
CTEMP=CW*CX(I+L)
CX(I+L)=CX(I)-CTEMP
CX(I)=CX(I)+CTEMP
END DO
END DO
L=ISTEP
END DO
IF(SIGNI.LE.0.0) THEN
SC=1./FLOAT(LX)
DO I=1,LX
CX(I)=CX(I)*SC
END DO
ENDIF
RETURN
END
SUBROUTINE DERIV(NTH,NWAY)
COMPLEX C,CON,CON2,CON3
COMMON/BLK/C(256,128),NX,NY,DELTA
C
C DERIV CALCULATES THE VALUES OF 2 QUADRANTS FOR THE NTH
C DERIVATIVE OF A (NX,NY) MATRIX. THESE VALUES ARE TO BE
C MUTLIPLIED BY THE WAVENUMBER SPECTRUM OF THE GIVEN FIELD.
C
C FOR A 1 DIMENSIONAL ARRAY SET NY=1
C
C NX =NUMBER OF COLUMNS IN THE DATA SET
C NY =NUMBER OF ROWS IN THE DATA SET
C NTH =THE ORDER OF DERIVATIVE DESIRED
C NWAY =THE DIRECTION THE DERIVATIVE IS TO BE TAKEN
C 1 VERTICAL DIRECTION
C 2 HORIZONTAL X DIRECTION
C 3 HORIZONTAL Y DIRECTION
C DELTA=GRID INTERVAL IN LENGTH UNITS
C
IF(NTH.LT.1.OR.NWAY.LT.1.OR.NWAY.GT.4) THEN
TYPE *,' INVALID PARAMETER FOR DERIV (NTH OR NWAY)'
STOP
ENDIF
NXX=NX+2
NX2=NX/2+1
NYY=NY+2
NY2=NY/2+1
RNXDEL=1.0/(FLOAT(NX)*DELTA)
RNYDEL=1.0/(FLOAT(NY)*DELTA)
C----------------------------------------------------------------
C TAKE VERTICAL DERIVATIVE IN WAVENUMBER DOMAIN
C----------------------------------------------------------------
IF(NWAY.EQ.1) THEN
TYPE 50,NTH,'VERTICAL '
50 FORMAT(' IDEAL FILTER, ORDER:',I2,3X,A12)
CON=(6.2831853,0.0)
C(1,1)=(0.0,0.0)
DO IX=2,NX2
FX2=FLOAT(IX-1)*RNXDEL
C(IX,1)=(CON*FX2)**NTH
C(NXX-IX,1)=C(IX,1)
END DO
IF(NY.GT.1) THEN
DO IY=2,NY2
IYY=NYY-IY
FY2=FLOAT(IY-1)*RNYDEL
CON2=(CON*FY2)**NTH
C(1,IY)=CON2
C(1,IYY)=CON2
FY2=FY2*FY2
DO IX=2,NX2
FX2=(FLOAT(IX-1)*RNXDEL)**2
CON2=(CON*SQRT(FX2+FY2))**NTH
C(NXX-IX,IY)=CON2
C(NXX-IX,IYY)=CON2
C(IX,IY)=CON2
C(IX,IYY)=CON2
END DO
END DO
ENDIF
C----------------------------------------------------------------
C TAKE HORIZONTAL Y DERIVATIVE IN WAVENUMBER DOMAIN
C----------------------------------------------------------------
ELSEIF(NWAY.EQ.3) THEN
TYPE 50, NTH,'HORIZONTAL Y'
CON=(0.0,6.2831853)*RNYDEL
DO IX=1,NX
C(IX,1)=(0.0,0.0)
END DO
IF(NY.EQ.1) THEN
TYPE *,'ERROR - CAN NOT CALCULATE ON 1 DIMENSION ARRAY'
ELSE
DO IY=2,NY2
IYY=NYY-IY
CON2=(FLOAT(IY-1)*CON)**NTH
CON3=CONJG(CON2)
DO IX=1,NX
C(IX,IYY)=CON3
C(IX,IY)=CON2
END DO
END DO
ENDIF
C----------------------------------------------------------------
C TAKE HORIZONTAL X DERIVATIVE IN WAVENUMBER DOMAIN
C----------------------------------------------------------------
ELSE
TYPE 50,NTH,'HORIZONTAL X'
CON=(0.0,6.2831853)*RNXDEL
DO IY=1,NY
C(1,IYP)=(0.0,0.0)
END DO
DO IX=2,NX2
CON2=(FLOAT(IX-1)*CON)**NTH
DO IY=1,NY
C(NXX-IX,IY)=CONJG(CON2)
C(IX,IY)=CON2
END DO
END DO
ENDIF
RETURN
END
SUBROUTINE CONTIN(ZCON)
COMPLEX H,CTEMP
COMMON/BLK/H(256,128),NX,NY,DELTA
C
C CONTIN COMPUTES TWO QUADRANTS (NX/2+1 BY NY) OF AN IDEAL
C UPWARD OR DOWNWARD CONTINUATION FILTER DIMENSIONED BY NX BY NY
C
C DELTA = GRID INTERVAL IN LENGTH UNITS
C ZCON = DEPTH OR HIGHT AT WHICH CONTINUATION IS DESIRED
C IN SAME UNITS AS DELTA
C IF ZCON IS NEGATIVE, FILTER WILL BE UPWARD CONTINUATION
C IF ZCON IS POSITIVE, FILTER WILL BE DOWNWARD CONTINUATION
C NX,NY = NUMBER OF COL,ROWS, MUST BE POWER OF 2
C
TYPE *,' CONTINUATION FILTER'
PI2Z=3.14159265*ZCON/DELTA
NXX=NX+2
NYY=NY+2
NX2=NX/2+1
RNX2=1.0/FLOAT(NX2)
NY2=NY/2+1
RNY2=1.0/FLOAT(NY2)
DO IY=1,NY2
IYY=NYY-IY
CON1=(FLOAT(IY-1)*RNY2)**2
X1=EXP(SQRT(CON1)*PI2Z)
H(1,IY)=CMPLX(X1,0.0)
IF(IY.NE.1) H(1,IYY)=H(1,IY)
DO IX=2,NX2
S=SQRT((FLOAT(IX-1)*RNX2)**2+CON1)
X1=EXP(S*PI2Z)
IXX=NXX-IX
CTEMP=CMPLX(X1,0.0)
H(IXX,IY)=CTEMP
H(IX,IY)=CTEMP
IF(IY.NE.1) THEN
H(IXX,IYY)=CTEMP
H(IX,IYY)=CTEMP
ENDIF
END DO
END DO
RETURN
END
SUBROUTINE BNDPAS(CCUTLO,CCUTHI,NPASS)
COMPLEX H,ZERO,ONE
COMMON/BLK/H(256,128),NX,NY,DELTA
C
C BNDPAS CALCULATES TWO QUADRANTS OF THE WAVENUMBER RESPONSE
C OF AN IDEAL BAND-PASS FILTER OF A (NX,NY) MATRIX.
C
C CCUTLO, CCUTHI RANGE THAT WILL EITHER BE PASSED OR REJECTED
C NPASS -1 TO REJECT, 1 TO PASS
C DELTA DATA GRID INTERVAL
C NX,NY NUMBER OF ROWS, COLUMNS
C
TYPE *,' BAND-PASS/REJECT FILTER'
WAVLEN=2.0*DELTA
FNQ1=1.0/WAVLEN
CUTHI=CCUTHI
IF(CUTHI.GT.FNQ1) CUTHI=FNQ1
CUTLO=CCUTLO
IF(CUTLO.LT.0.0) CUTLO=0.0
NXX=NX+2
NYY=NY+2
NX2=(NX/2)+1
NY2=(NY/2)+1
IF(NPASS.EQ.1) THEN
ZERO=(0.0,0.0)
ONE=(1.0,0.0)
ELSE
ZERO=(1.0,0.0)
ONE=(0.0,0.0)
ENDIF
RHIY2=(FLOAT(NY2)*CUTHI*WAVLEN)**2
RLOWY2=(FLOAT(NY2)*CUTLO*WAVLEN)**2
CLOWX=FLOAT(NX2)*WAVLEN*CUTLO
CHIX=FLOAT(NX2)*WAVLEN*CUTHI
C
C SET ROW 1 AND COLUMN 1 TO 0.0
C
DO IY=1,NY
H(1,IY)=(1.,0.)
END DO
DO IX=1,NX
H(IX,1)=(1.,0.)
END DO
C----------------------------------------------------------------
C SET SELECT RANGE TO "ONE", REST TO "ZERO"
C----------------------------------------------------------------
C
C X**2 Y**2 A = CLOWX OR CHIX
C FOR ---- + ---- = 1 B**2 = RLOWY2 OR RHIY2
C A**2 B**2
C
C OR X = A * SQRT(1-Y**2/B**2)
C
DO IY=2,NY2
Y2=FLOAT(IY-1)**2
IYY=NYY-IY
IF(Y2.LT.RLOWY2) THEN
MINX=CLOWX*SQRT(1.0-Y2/RLOWY2)+1.0001
ELSE
MINX=0.0
ENDIF
IF(Y2.LE.RHIY2) THEN
MAXX=CHIX*SQRT(1.0-Y2/RHIY2)+1.0001
ELSE
MAXX=0.0
ENDIF
C
IF(MINX.GE.2) THEN
DO IX=2,MINX-1
H(IX,IY)=ZERO
H(IX,IYY)=ZERO
H(NXX-IX,IYY)=ZERO
H(NXX-IX,IY)=ZERO
END DO
ENDIF
C
IF(MAXX.GE.2) THEN
IF(MINX.LT.2) MINX=2
DO IX=MINX,MAXX
H(IX,IY)=ONE
H(IX,IYY)=ONE
H(NXX-IX,IYY)=ONE
H(NXX-IX,IY)=ONE
END DO
ELSE
MAXX=1
ENDIF
IF(MAXX.LT.NX2) THEN
DO IX=MAXX+1,NX2
H(IX,IY)=ZERO
H(IX,IYY)=ZERO
H(NXX-IX,IYY)=ZERO
H(NXX-IX,IY)=ZERO
END DO
ENDIF
END DO
RETURN
END
SUBROUTINE STRIKE(AANG1,AANG2,NNPASS)
COMPLEX H,ZERO,ONE
COMMON/BLK/H(256,128),NX,NY,DELTA
DATA D2RAD,DEG90/0.0174532,1.5707963/
C
C STRIKE CREATES A STRIKE SENSITIVE FILTER (FAN FILTER)
C FOR 2 QUADRANTS OF THE (NX,NY) MATRIX
C
C ANGLES ARE MEASURED IN DEGREES CLOCKWISE FROM NORTH
C
C ANG1 SMALLEST ANGLE GE 0.0 AND LT ANG2
C ANG2 LARGEST ANGLE GT ANG1 AND LE 180.
C NPASS STATES IF DATA IS PASSED OR REJECTED BETWEEN ANGLES
C -1 IF REJECT AZIMUTHS BETWEEN ANGLES
C 1 IF PASS AZIMUTHS BETWEEN ANGLES
C NX NY NUMBERS OF ROWS AND COLUMNS
C
TYPE *,' STRIKE PASS/REJECT FILTER'
ANG1=AANG1
ANG2=AANG2
NPASS=NNPASS
NX2=NX/2+1
NY2=NY/2+1
NXX=NX+2
NYY=NY+2
XY=FLOAT(NX2)/FLOAT(NY2)
IF(NPASS.NE.-1) THEN
ZERO=(0.0,0.0)
ONE=(1.0,0.0)
ELSE
ZERO=(1.0,0.0)
ONE=(0.0,0.0)
ENDIF
C----------------------------------------------------------------
C SET MATRIX TO "ZERO"
C----------------------------------------------------------------
DO IY=1,NY
DO IX=1,NX
H(IX,IY)=ZERO
END DO
END DO
C----------------------------------------------------------------
C COMPUTE PARAMETERS FOR SOUTHWEST QUADRANT OF MATRIX
C----------------------------------------------------------------
IF(ANG2.LT.90.0) THEN
IYMAX=0
ELSE
IF(ANG1.LE.90.0) THEN
A1=0.0
TA1=0.0
IYMAX=NY2
ADA1=1.5
ELSE
A1=ATAN(TAN((ANG1-90.0)*D2RAD)*XY)
TA1=TAN(A1)
IYMAX=FLOAT(NX2)/TA1+1.0
IF(IYMAX.GT.NY2) IYMAX=NY2
ADA1=1.0
ENDIF
IF(ANG2.GE.180.0) THEN
A2=DEG90
TA2=0.0
ADA2=FLOAT(NX2)+0.5
ELSE
A2=ATAN(TAN((ANG2-90.0)*D2RAD)*XY)
TA2=TAN(A2)
ADA2=1.0
ENDIF
ENDIF
IF(ANG1.GT.90.0) THEN
IYMAXX=0
ELSE
C----------------------------------------------------------------
C COMPUTE PARAMETERS FOR SOUTHEAST QUADRANT OF MATRIX
C----------------------------------------------------------------
IF(ANG1.LE.0.0) THEN
A1=0.0
TTA1=0.0
ADDA1=FLOAT(NX2)-2.5
ELSE
IF(ANG1.GT.90.0) THEN
A1=0.0
A2=0.0
TTA1=0.0
TTA2=0.0
ADDA1=0.0
ADDA2=0.0
IYMAX=NY2
GOTO 100
ELSE
A1=ATAN(TAN(ANG1*D2RAD)/XY)
TTA1=1.0/TAN(A1)
ADDA1=0.0
ENDIF
ENDIF
IF(ANG2.GE.90.0) THEN
IYMAXX=NY2
TTA2=0.0
ADDA2=0.0
ELSE
A2=ATAN(TAN(ANG2*D2RAD)/XY)
TTA2=1.0/TAN(A2)
ADDA2=0.0
IYMAXX=ABS(FLOAT(NX2)/TTA2+1.5)
IF(IYMAXX.GT.NY2) IYMAXX=NY2
ENDIF
ENDIF
C================================================================
C CALCULATE THE FILTER COEFFICENTS
C================================================================
100 NYMAX=AMAX0(IYMAX,IYMAXX)
H(1,1)=(1.0,0.0)
DO IY=1,NY2
DO IX=1,NX
H(IX,IY)=ZERO
END DO
IF(IY.LE.NYMAX) THEN
Y=FLOAT(IY-1)
C----------------------------------------------------------------
C DEFINE SOUTHWEST QUADRANT
C----------------------------------------------------------------
IF(IYMAX.GE.IY) THEN
MIN=Y*TA1+ADA1
MAX=Y*TA2+ADA2
IF(MAX.GT.NX2) MAX=NX2
IF(MIN.LE.MAX) THEN
DO IX=MIN,MAX
H(IX,IY)=ONE
END DO
ENDIF
ENDIF
C----------------------------------------------------------------
C DEFINE SOUTHEAST QUADRANT
C----------------------------------------------------------------
IF(IYMAXX.GE.IY) THEN
MIN=NX-(Y*TTA2+ADDA2)
IF(MIN.GT.NX) MIN=NX
MAX=NX-(Y*TTA1+ADDA1)
IF(MAX.LT.NX2+1) MAX=NX2+1
DO IX=MAX,MIN
H(IX,IY)=ONE
END DO
ENDIF
ENDIF
C----------------------------------------------------------------
C USE ANTI SYMMETRY TO DEFINE QUADRANTS 2 AND 3
C----------------------------------------------------------------
IF(IY.NE.1.AND.IY.NE.NY2) THEN
IYY=NYY-IY
DO IX=2,NX2
H(IX,IYY)=H(NXX-IX,IY)
H(NXX-IX,IYY)=H(IX,IY)
END DO
ENDIF
END DO
RETURN
END
54 69 6D 20 46 6F 67 61 72 74 79 20 6C 6F 76 65 73 20
53 74 65 76 65 20 52 69 6F 72 64 61 6E 21
[Editors note: As it is bound and sitting on some shelf in the Geoscience Library, I placed the previous coded message in various locations in my thesis. As I was 25 and only recently out, this was making quite a statement. $10 or a big kiss, which every you prefer, to the first person who emails me the correct uncoded message.]
One of the easer methods to interactively run a program is through the use of menus. The entire screen of the video terminal can be used to list options and parameters of the program. The user can enter commands to select options, or move a cursor to select desired parameters. Menus can be set up where each input parameter is labeled and displayed. On VT-100 compatible terminals, the PF keys can be used for special functions. The arrow keys can be used to move the cursor on the screen to the desired parameter, and new data can be entered. The subroutines listed here allow easy access to the VAX/VMS Input/Output procedures in the Run Time Library, and with them, programs can be written in the menu format. For a more detailed description of the Input/ Output procedures, see chapter 3 of the VAX-11 Run Time Library User's Guide in VAX/VMS (version 3.2) volume 6A. It is understood that VAX/VMS version 4.0 will have a much more detailed set of screen input/output routines as part of its run time library.
Listed here are brief descriptions of each of the subroutines used to create and manipulate menus on VT-100 compatible screens. The following section lists the source code for these routines.
If the user wanted to print a real number in F10.3 format, LNUM would equal 10 and LDEC would equal 3.
SUBROUTINE VT_SETUP
INTEGER SYS$ASSIGN
INTEGER*2 IFLAG
CHARACTER*4 DEV
COMMON/IOCHAN/IO_CHAN,IO_EFN,IO_FUN
EXTERNAL IO$_READPROMPT,IO$M_TRMNOECHO,IO$M_ESCAPE
DATA DEV/'TT: '/
IDSW=SYS$ASSIGN(DEV,IO_CHAN,,)
CALL LIB$GET_EF(IO_EFN)
IF(IDSW.NE.1) CALL LIB$STOP(%VAL(IDSW))
ISTATUS=LIB$SCREEN_INFO(IFLAG)
IO_FUN=%LOC(IO$_READPROMPT)+
> %LOC(IO$M_TRMNOECHO)+%LOC(IO$M_ESCAPE)
RETURN
END
SUBROUTINE ERASE_VT
ISTATUS=LIB$ERASE_PAGE(1,1)
RETURN
END
SUBROUTINE MOVE_VT(I,J)
ISTATUS=LIB$SET_CURSOR(I,J)
RETURN
END
SUBROUTINE ERASE_LINE(I,J)
ISTATUS=LIB$ERASE_LINE(I,J)
RETURN
END
SUBROUTINE TEXT_VT(TEXT,I,J,IFLAG)
CHARACTER*(*) TEXT
LENGTH=LEN(TEXT)
ISTATUS=LIB$PUT_SCREEN(TEXT(1:LENGTH),I,J,IFLAG)
RETURN
END
SUBROUTINE PUTRNUM(RNUM,I,J,IFLAG,M,N)
CHARACTER*20 TEXT
WRITE(TEXT,'(6X,F14.4)') RNUM
IE=16+N
IF(IE.GT.20) IE=20
IB=IE-M
ISTATUS=LIB$PUT_SCREEN(TEXT(IB:IE),I,J,IFLAG)
RETURN
END
SUBROUTINE PUTINUM(INUM,I,J,IFLAG,M)
CHARACTER*11 TEXT
WRITE(TEXT,'(2X,I8)') INUM
K=10-M
ISTATUS=LIB$PUT_SCREEN(TEXT(K:10),I,J,IFLAG)
RETURN
END
SUBROUTINE ASKQ(PROMPT,ANSWER,I,J)
CHARACTER*(*) PROMPT,ANSWER
ISTATUS=LIB$ERASE_LINE(I,J)
ISTATUS=LIB$SET_CURSOR(I,J)
ISTATUS=LIB$GET_SCREEN(ANSWER,PROMPT,LENGTH)
ISTATUS=STR$UPCASE(ANSWER,ANSWER)
RETURN
END
SUBROUTINE ASKQ2(ANSWER,I,J)
CHARACTER*(*) ANSWER
ISTATUS=LIB$SET_CURSOR(I,J)
ISTATUS=LIB$GET_SCREEN(ANSWER)
ISTATUS=STR$UPCASE(ANSWER,ANSWER)
RETURN
END
SUBROUTINE ASKQPF(PROM,ANSWER,IP,JP,IPFVAL)
INTEGER SYS$QIOW
INTEGER*2 IOSB(4)
CHARACTER PROM*(*),ANSWER*(*),TEXT*32,CHAR2,BRACKET,OCHAR
BYTE BUFPROM(80),BUF(80)
COMMON/IOCHAN/IO_CHAN,IO_EFN,IO_FUNC
DATA IA,IP,OCHAR/65,80,'O'/ ! ASCII FOR A, P, O
BRACKET=CHAR(91) ! ASCII CODE FOR
C ! LEFT SQAURE BRACKET
LEN_PROM=LEN(PROM)
ISTATUS=LIB$ERASE_LINE(IP,JP)
ISTATUS=LIB$SET_CURSOR(IP,JP)
DO I=1,LEN_PROM
BUFPROM(I)=ICHAR(PROM(I:I))
END DO
IDSW1=SYS$QIOW(%VAL(IO_EFN),%VAL(IO_CHAN),%VAL(IO_FUNC),
> IOSB,,,BUF,%VAL(80),,,BUFPROM,%VAL(LEN_PROM))
IF(IDSW1.GT.1) CALL LIB$STOP(%VAL(IDSW1))
LENOUT=IOSB(2)
ANSWER=CHAR(0)
DO I=1,LENOUT
ANSWER(I:I)=CHAR(BUF(I))
END DO
IPFVAL=0
IF(BUF(LENOUT+1).EQ.27) THEN
CHAR2=CHAR(BUF(LENOUT+2))
I=BUF(LENOUT+3)
IF(CHAR2.EQ.BRACKET) IPFVAL=I-IA+5
IF(CHAR2.EQ.OCHAR) IPFVAL=I-IP+1
ENDIF
ISTATUS=STR$UPCASE(ANSWER,ANSWER)
RETURN
END
SUBROUTINE TXT2INT(TEXT,INUM)
CHARACTER*(*) TEXT
L=INDEX(TEXT,'.')
IF(L.NE.0) THEN
LENGTH=LEN(TEXT)
DO LI=L,LENGTH
TEXT(LI:LI)=' '
END DO
ENDIF
READ(TEXT,'(BNI10)') INUM
RETURN
END
SUBROUTINE TXT2INT2(TEXT,INUM)
INTEGER*2 INUM
CHARACTER*(*) TEXT
L=INDEX(TEXT,'.')
IF(L.NE.0) THEN
LENGTH=LEN(TEXT)
DO LI=L,LENGTH
TEXT(LI:LI)=' '
END DO
ENDIF
READ(TEXT,'(BNI5)') INUM
RETURN
END
SUBROUTINE TXT2R(TEXT,RNUM)
CHARACTER*(*) TEXT
L=INDEX(TEXT,' ')
IF(INDEX(TEXT,'.').EQ.0) TEXT(L:L+1)='. '
READ(TEXT,'(BNF10.4)') RNUM
RETURN
END
SUBROUTINE BELL
CHARACTER*1 RING
DATA I,J,RING/25,1,7/
ISTATUS=LIB$PUT_SCREEN(RING,I,J)
RETURN
END
LOGICAL*2 FUNCTION LNUMBER(TEXT)
INTEGER F
LOGICAL*2 FIRST
CHARACTER*(*) TEXT
LENGTH=LEN(TEXT)
NUMBER=.FALSE.
FIRST=.TRUE.
DO I=1,LENGTH
F=ICHAR(TEXT(I:I))
IF(F.NE.32) THEN
IF(FIRST) THEN
NUMBER=F.EQ.43.OR.(F.GE.45.AND.F.LE.57.AND.F.NE.47)
FIRST=.FALSE.
ELSE
NUMBER=F.GE.46.AND.F.LE.57.AND.F.NE.47
ENDIF
IF(.NOT.NUMBER) RETURN
ENDIF
END DO
RETURN
END
SUBROUTINE SCROLL(I,J)
COMMON/IOCHAN/IO_CHAN,IO_EFN,IO_FUN,ISCROLL,JSCROLL
ISCROLL=I
JSCROLL=J
ISTATUS=LIB$SET_SCROLL(ISCROLL,JSCROLL)
RETURN
END
SUBROUTINE UP_VT
COMMON/IOCHAN/IO_CHAN,IO_EFN,IO_FUN,ISCROLL,JSCROLL
ISTATUS=LIB$SET_CURSOR(JSCROLL,1)
ISTATUS=LIB$UP_SCROLL()
RETURN
END
SUBROUTINE DOWN_VT
COMMON/IOCHAN/IO_CHAN,IO_EFN,IO_FUN,ISCROLL,JSCROLL
ISTATUS=LIB$SET_CURSOR(ISCROLL,1)
ISTATUS=LIB$DOWN_SCROLL()
RETURN
END
54 69 6D 20 46 6F 67 61 72 74 79 20 6C 6F 76 65 73 20 53 74 65 76 65 20 52 69 6F 72 64 61 6E 21
Editors note: As it is bounded and sitting on some shelf in the Geoscience Library at USC, I placed the above coded message in various locations in my thesis. As I was 25 and only recently out, this was making quite a statement. Can you figure it out ?]
Next Chapter: Appendix A, section 5
Previous Chapter: Bibliography