CTITLESARAMRH - 2D DECON FOR SPACE FREQUENCY DOMAIN RADON TRANSFORM 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA 00020000 CA AUTHOR CHUCK MOSHER 00030000 CA DESIGNER CHUCK MOSHER 00040000 CA LANGUAGE FORTRAN 77 00050000 CA SYSTEM IBM AND CRAY 00060000 CA WRITTEN 09-01-87 00070000 C REVISED 03-06-89 JJC - FOR SPARC AND CRAY PRODUCTION. 00080000 C REVISED 04-20-90 ESN - FOR CRAY CONVERSION. 00090000 C REVISED 02-13-92 ESN - COPIED FROM SARAMRD FOR 00100000 C PERFORMANCE IMPROVEMENTS. 00110000 CA 00120000 CA CALL SARAMRH(NF,SW,NP,DP,DX,RA,NTRC,D,WNOIS, 00130000 CA + R,C,H,A,B,IPR,KPBUGF) 00140000 CA 00150000 CA PURPOSE: 2-D DECON FOR SPACE-FREQUENCY DOMAIN RADON TRANSFORMS 00160000 CA 00170000 CA INPUT NF NUMBER OF FREQUENCIES TO BE PROCESSED I400180000 CA INPUT SW FIRST FREQUENCY SAMPLE R400190000 CA INPUT NP NUMBER OF P'S I400200000 CA INPUT DP P INCREMENT R400210000 CA INPUT DX MAXIMUM TRACE SEPARATION R400220000 CA I/O RA INPUT RADON TRANSFORM, OUTPUT DECONVOLVED TRANSFORMC800230000 CA INPUT NTRC NUMBER OF TRACES I400240000 CA INPUT D TRACE DISTANCE ARRAY R400250000 CA INPUT WNOIS WHITE NOISE LEVEL FOR DECON R400260000 CA INPUT R COMPLEX WORK VECTORS, LENGTH NTRC C800270000 CA INPUT C COMPLEX WORK VECTORS, LENGTH NTRC C800280000 CA INPUT H COMPLEX WORK VECTORS, LENGTH NTRC C800290000 CA INPUT A COMPLEX WORK VECTORS, LENGTH NTRC C800300000 CA INPUT B COMPLEX WORK VECTORS, LENGTH NTRC C800310000 CA INPUT IPR PRINT UNIT NUMBER, 0 FOR NO PRINT I400320000 C--------------------------------------------------------------------- 00330000 C 00340000 SUBROUTINE SARAMRH(NF,SW,NP,DP,DX,RA,NTRC,D,WNOIS, 00350000 + R,C,H,A,B,IPR,KPBUGF) 00360000 C 00370000 IMPLICIT INTEGER (A-Z) 00380000 C 00390000 COMPLEX RA(NF,*),R(1),C(1),H(1),A(1),B(1),RSUM,RDC 00400000 COMPLEX CZERO 00410000 DIMENSION FREQS(10),D(1) 00420000 C 00430000 REAL D 00440000 REAL DP 00450000 REAL DX 00460000 REAL FACT 00470000 REAL FREQS 00480000 REAL PI2 00490000 REAL RHZ 00500000 REAL RPOW 00510000 REAL SW 00520000 REAL W 00530000 REAL WNOIS 00540000 REAL WP 00550000 REAL XNOIS 00560000 C 00570000 DATA PI2/6.2831853/,CZERO /(0.,0.)/ 00580000 C 00590000 DATA IFIRST / 0 / 00600000 INTEGER SCINDX (8000) 00610000 REAL ACINDX (8000) 00620000 REAL SCTAB (50001) 00630000 C 00640000 COMMON /SYSTEM/ SYSTEM 00650000 C 00660000 IF (IFIRST .EQ. 0) THEN 00670000 IFIRST = 1 00680000 C -- NOTE: NSCPTS MUST BE DIVISIBLE BY 4 (FOR COS LOOKUP) 00690000 NSCPTS = 40000 00700000 C 00710000 INDSIN = 1 00720000 DO 1 I = 1, NSCPTS 00730000 WP = ((I-1)*PI2) / NSCPTS 00740000 SCTAB(INDSIN+I-1) = SIN(WP) 00750000 1 CONTINUE 00760000 C 00770000 INDCOS = INDSIN + NSCPTS/4 00780000 DO 2 I = 1, NSCPTS/4+1 00790000 SCTAB(INDSIN+NSCPTS+I-1) = SCTAB(INDSIN+I-1) 00800000 2 CONTINUE 00810000 ENDIF 00820000 IF (NTRC .GT. 8000) CALL XDUMPX (NTRC) 00830000 IF (NTRC .LT. 1) CALL XDUMPX (NTRC) 00840000 C 00850000 RHZ = 1./PI2 00860000 FACT = PI2/(DX*DP) 00870000 IF (KPBUGF.NE.0) 00880000 + WRITE (IPR,'(/,'' P DECON: FREQUENCIES PROCESSED '')') 00890000 C 00900000 C CONSTRUCT ARRAYS FOR DECON 00910000 C ARRAY USAGE: 00920000 C R = FIRST ROW OF TOEPPLITZ MATRIX 00930000 C C = TEMPORARY STORAGE OF ONE COLUMN OF DATA ARRAY 00940000 C H = OUTPUT RESULT FROM LEVINSON RECURSION 00950000 C A,B = WORK ARRAYS FOR RECURSION 00960000 C 00970000 C OUTER LOOP IS ON FREQUENCY 00980000 C 00990000 DO 40 I=2,NF 01000000 C 01010000 C SET FREQUENCY 01020000 C 01030000 W = SW*(I-1) 01040000 C 01050000 C COMPUTE NUMBER OF UN-ALIASED P VALUES 01060000 C 01070000 NPA = 1 + IFIX(FACT/W) 01080000 IF (NPA.GT.NP) NPA = NP 01090000 C 01100000 C INITIALIZE DC LEVEL AND TOTAL POWER FOR NOISE FACTOR COMPUTATION 01110000 C 01120000 RDC = (0.,0.) 01130000 RPOW = 0. 01140000 C 01150000 C LOOP OVER NUMBER OF P VALUES 01160000 C 01170000 DO 20 K=1,NPA 01180000 RSUM = (0.,0.) 01190000 C 01200000 C P TIMES FREQUENCY - NOTE THAT IT IS POSITIVE 01210000 C 01220000 WP = W*DP*(K-1) 01230000 C 01240000 C LOOP OVER TRACES AND CONSTRUCT ELEMENT OF DCON OPERATOR 01250000 C 01260000 C DO 10 J=1,NTRC 01270000 C RSUM = RSUM + CMPLX(COS(WP*D(J)),SIN(WP*D(J))) 01280000 C 10 CONTINUE 01290000 DO 10 J=1,NTRC 01300000 ACINDX(J) = (WP*D(J)) / PI2 01310000 10 CONTINUE 01320000 DO 12 J=1,NTRC 01330000 ACINDX(J) = ACINDX(J) - INT(ACINDX(J)) 01340000 12 CONTINUE 01350000 DO 13 J=1,NTRC 01360000 IF (ACINDX(J) .LT. 0.0) ACINDX(J) = ACINDX(J) + 1.0 01370000 SCINDX(J) = ACINDX(J)*NSCPTS + 0.5 01380000 13 CONTINUE 01390000 C DO 14 J=1,NTRC 01400000 C IF (SCINDX(J) .LT. 0) CALL XDUMPX (J,SCINDX(J)) 01410000 C IF (SCINDX(J) .GT. 40000) CALL XDUMPX (J,SCINDX(J)) 01420000 C 14 CONTINUE 01430000 IF (S1CPCH(SYSTEM,1,'CRAY',1,4) .EQ. 0) THEN 01440000 DO 15 J=1,NTRC 01450000 RSUM = RSUM + CMPLX(SCTAB(INDCOS+SCINDX(J)), 01460000 * SCTAB(INDSIN+SCINDX(J))) 01470000 15 CONTINUE 01480000 ELSE 01490000 J = 0 01500000 16 J = J + 1 01510000 IF (J .GT. NTRC) GO TO 17 01520000 IF (J .GT. 8000) CALL XDUMPX 01530000 RSUM = RSUM + CMPLX(SCTAB(INDCOS+SCINDX(J)), 01540000 * SCTAB(INDSIN+SCINDX(J))) 01550000 GO TO 16 01560000 17 CONTINUE 01570000 ENDIF 01580000 C 01590000 C SUM DC AND TOTAL POWER OF OPERATOR 01600000 C 01610000 R(K) = RSUM 01620000 RPOW = RPOW + REAL(RSUM)**2 + AIMAG(RSUM)**2 01630000 RDC = RDC + RSUM 01640000 C 01650000 C MOVE COLUMN OUT OF DATA ARRAY 01660000 C 01670000 C(K) = RA(I,K) 01680000 20 CONTINUE 01690000 C 01700000 C ADD NOISE TO FIRST VALUE OF R 01710000 C 01720000 IF (RPOW .NE. 0.0) THEN 01730000 XNOIS = 1. + .01*WNOIS*REAL(RDC*CONJG(RDC))/RPOW 01740000 ELSE 01750000 XNOIS = 1.0 01760000 ENDIF 01770000 IF (R(1) .NE. CZERO) THEN 01780000 R(1) = XNOIS*R(1) 01790000 ELSE 01800000 R(1) = CMPLX(.01*WNOIS,.01*WNOIS) 01810000 ENDIF 01820000 C 01830000 C SOLVE SET OF EQUATIONS C = R*H FOR H 01840000 C 01850000 CALL ARSET(H,2*NP,0) 01860000 CALL SACLEVN(NPA,R,C,H,A,B) 01870000 C 01880000 C MOVE SOLUTION AT THIS FREQUENCY BACK INTO RA 01890000 C 01900000 DO 30 J=1,NP 01910000 RA(I,J) = H(J) 01920000 30 CONTINUE 01930000 C 01940000 C PRINT FREQUENCIES PROCESSED 01950000 C 01960000 IF (KPBUGF.NE.0) THEN 01970000 L = MOD(I-2,10) + 1 01980000 FREQS(L) = W*RHZ 01990000 IF (L.EQ.10) WRITE (IPR,9000) FREQS 02000000 9000 FORMAT (' ','PROCESSED FREQUENCIES:',/,10F7.1) 02010000 ENDIF 02020000 40 CONTINUE 02030000 IF (KPBUGF.NE.0 .AND. L.NE.10) WRITE (IPR,9000) 02040000 + (FREQS(I),I=1,L) 02050000 RETURN 02060000 END 02070000