CTITLESAMGFK -- FK MIGRATION WITH STRETCH USING WIENER FILTERS 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR C. J. SICKING 00000200 CA DESIGNER C. J. SICKING 00000300 CA LANGUAGE FORTRAN 00000400 CA SYSTEM IBM AND CRAY 00000500 CA WRITTEN SEPTEMBER, 1980 00000600 C REVISED 03-27-81 DJP. CHANGED NAME FROM SAM2FK, ADDED DIP 00000700 C FILTERING AND 3-D PROCESSING CAPABILITY. 00000800 C REVISED 02-18-83 ESN. ADDED PRE-STACK CAPABILITY. 00000900 C REVISED 12-10-84 LBL. MADE THIS CODE RAN ON THE CRAY. 00001000 C REVISED 10-30-85 ESN. CHANGE '&' TO '*' IN ALTERNATE RETURNS.00001100 C REVISED 11-13-89 RDK. FOR CRAY CFT77 COMPATIBILITY. 00001202 CA 00001300 CA 00001400 CA CALL SAMGFK (NX, NT, NXP, NTP, NWP, NEXPT, NEXPX, DX, DT, ND, DD, 00001500 CA DTOT, TTOT, VC, IOADR, IOAD2, THL, NXBUF, WAV, C1, 00001600 CA Y, Z, WBUF, COEF, FRACN, KPPRNT, DIPPS, DIPRJ, NXM, 00001700 CA MODE, ILINE, IDISK, JF, CDPFLD, IER) 00001800 CA 00001900 CA IN/OUT ARGUMENT TYPE DESCRIPTION 00002000 CA 00002100 CA IN NX I4 NUMBER OF SAMPLES IN SPACE 00002200 CA IN NT I4 NUMBER OF SAMPLES IN TIME 00002300 CA IN NXP I4 NUMBER OF POINTS IN SPACE FFT 00002400 CA IN NTP I4 NUMBER OF POINTS IN TIME FFT 00002500 CA IN NWP I4 NTP/2 +1 00002600 CA IN NEXPT I4 POWER OF 2 FOR TIME FFT 00002700 CA IN NEXPX I4 POWER OF 2 FOR SPACE FFT 00002800 CA IN DX R4 DELTA X 00002900 CA IN DT R4 DELTA T 00003000 CA IN ND I4 NUMBER OF POINTS IN RESAMPLED TRACE 00003100 CA IN DD R4 SAMPLING RATE OF RESAMPLED TRACES 00003200 CA IN DTOT R4 DATA ARRAY FOR DEPTH TO TIME RESAMPLING 00003300 CA IN TTOD R4 DATA ARRAY FOR TIME TO DEPTH RESAMPLING 00003400 CA IN VC R4 RMS VELOCITY VALUE AT SAMPLE NT-1 00003500 CA IN IOADR I4 ID FOR WORKFILE #1 00003600 CA IN IOAD2 I4 ID FOR WORKFILE #2 00003700 CA IN THL R4 LENGTH OF HEADER 00003800 CA IN NXBUF I4 LENGTH OF 2D MEMORY BUFFER IN SPACE 00003900 CA IN WAV R4 2D COMPLEX ARRAY - STORAGE AREA FOR X-T 00004000 CA SECTION AND F-K SPECTRUM 00004100 CA IN C1 R4 COMPLEX WORK AREA 00004200 CA IN Y R4 REAL WORK AREA 00004300 CA IN Z R4 REAL WORK AREA 00004400 CA IN WBUF R4 ARRAY OF LENGTH NT+THL 00004500 CA IN COEF R4 WORK ARRAY FOR 99 INTERPOLATION FILTERS 00004600 CA IN FRACN R4 MAXIMUM FREQUENCY AS A FRACTION OF THE 00004700 CA NYQUIST FREQUENCY 00004800 CA IN KPPRNT I4 LOGICAL UNIT NUMBER FOR PRINTER OUTPUT 00004900 CA IN DIPPS I4 PASS DIP ANGLE IN DEGREES 00005000 CA IN DIPRJ I4 REJECT DIP ANGLE IN DEGREES 00005100 CA IN NXM I4 MAXIMUM NUMBER OF CDP ON ANY 3-D LINE 00005200 CA IN MODE I4 INDICATOR OF MIGRATION DIRECTION 00005300 CA 1 = X DIRECTION 00005400 CA 2 = Y DIRECTION 00005500 CA 3 = BOTH DIRECTIONS, X DIR. CURRENTLY 00005600 CA 4 = BOTH DIRECTIONS, Y DIR. CURRENTLY 00005700 CA IN ILINE I4 SEQUENTIAL LINE NUMBER 00005800 CA IN IDISK I4 WORKFILE USEAGE FLAG : 0 = USE WORKFILE 00005900 CA IN JF I4 CDP PRE-STACK TRACE NUMBER 00006000 CA IN CDPFLD I4 CDP PRE-STACK MAXIMUM FOLD 00006100 CA OUT IER I4 FOURIER TRANSFORM ERROR FLAG 00006200 CA 00006300 CA 00006400 CA PERFORMS FK MIGRATION WITH STRETCH USING WIENER FILTERS IN 00006500 CA ALL INTERPOLATIONS. 00006600 CA 00006700 CA 00006800 C 00006900 SUBROUTINE SAMGFK (NX, NT, NXP, NTP, NWP, NEXPT, NEXPX, DX, DT, 00007000 * ND, DD, DTOT, TTOD, VC, IOADR, IOAD2, THL, 00007100 * NXBUF, WAV, C1, Y, Z, WBUF, COEF, FRACN, 00007200 * KPPRNT, DIPPS, DIPRJ, NXM, MODE, ILINE, 00007300 * IDISK, JF, CDPFLD, IER) 00007400 C 00007500 COMPLEX C1 ( 1) 00007600 COMPLEX WAV (NXBUF, 1) 00007700 C 00007800 COMPLEX WA 00007900 C 00008000 REAL A ( 14) 00008100 REAL COEF (12, 99) 00008200 REAL DTOT ( 1) 00008300 REAL FF ( 13) 00008400 REAL G ( 12) 00008500 REAL RW ( 12) 00008600 REAL TTOD ( 1) 00008700 REAL WBUF ( 1) 00008800 REAL Y ( 1) 00008900 REAL Z ( 1) 00009000 C 00009100 INTEGER CDPFLD 00009200 INTEGER DIPPS 00009300 INTEGER DIPRJ 00009400 INTEGER THL 00009500 C 00009600 C COMPUTE PARAMETERS NEEDED FOR THE F-K MAPPING STEP 00009700 C 00009800 IER = 0 00009900 PI = 3.14159265 00010000 DKD = 2.0 * PI / (NTP * DD) 00010100 DKZ = 2.0 * DKD / VC 00010200 DKX = 2.0 * PI / (NXP * DX) 00010300 C 00010400 C CHANGE DIP FILTER PARAMETERS FROM DEGREES TO RADIANS 00010500 C 00010600 DIPW = DIPRJ 00010700 DIPS = DIPPS 00010800 DIPS = PI * DIPS / 180. 00010900 DIPW = PI * DIPW / 180. 00011000 C 00011100 C FORWARD 2D FFT ***************************************************00011200 C 00011300 C ND1 TO NTP WILL BE ZEROED IN THE TRACE ARRAY BEFORE THE TIME FFT 00011400 C NX1 TO NXP WILL BE ZEROED IN THE SPACE ARRAY BEFORE THE SPACE FFT 00011500 C 00011600 ND1 = ND + 1 00011700 NX1 = NX + 1 00011800 C 00011900 C GENERATE THE INTERPOLATION COEFFICIENTS FOR STRETCHING THE TRACES.00012000 C THE COEFFICIENTS ARE COMPUTED HERE FOR USE LATER. 00012100 C 00012200 LR = 10 00012300 IF (FRACN .GT. 0.62) LR = 12 00012400 BW = FRACN 00012500 IF (FRACN .GT. 0.85) BW = 0.85 00012600 CALL MCOFGN (LR, BW, RW, G, FF, A, COEF) 00012700 C 00012800 C FOURIER TRANSFORM, T --> W FOR ALL TRACES 00012900 C 00013000 C THE TRACES ARE READ FROM THE TRACE SAVE FILE, RESAMPLED, FOURIER 00013100 C TRANSFORMED AND STORED IN WAV. AFTER NXBUF TRACES ARE READ 00013200 C WAV IS WRITTEN TO THE FFT SAVE FILE AND THE NEXT GROUP OF 00013300 C NXBUF TRACES ARE READ FROM THE TRACE SAVE FILE. IF IDISK IS 00013400 C EQUAL TO ONE, THE 2D FFT WILL FIT IN THE MEMORY BUFFER WAV 00013500 C AND THE DISC STORAGE IS BYPASSED. 00013600 C 00013700 IOUTR = 1 00013800 IKX = 0 00013900 KX = 1 00014000 NTX = NXBUF 00014100 IF (NX .GE. NXBUF) GO TO 10 00014200 NTX = NX 00014300 C 00014400 C THIS DO LOOP IS EXECUTED NXP / NXBUF TIMES 00014500 C 00014600 10 DO 60 IX = KX, NTX 00014700 IKX = IKX + 1 00014800 C 00014900 C GET INPUT TRACE 00015000 C 00015100 IF (MODE .EQ. 1) IOREC = ((ILINE-1) * CDPFLD + JF - 1) * NXM + IX 00015200 IF (MODE .EQ. 3) IOREC = ((ILINE-1) * CDPFLD + JF - 1) * NXM + IX 00015300 IF (MODE .EQ. 2) IOREC = ((IX-1) * CDPFLD + JF - 1) * NXM + ILINE 00015400 IF (MODE .EQ. 4) IOREC = ((IX-1) * CDPFLD + JF - 1) * NXM + ILINE 00015500 CALL FORDSD ( IOADR, IOREC, WBUF) 00015600 CALL ARMVE ( WBUF(THL+1), Z, NT) 00015700 IF (MODE .EQ. 4) CALL ARMVE (Z, Y, NT) 00015800 IF (MODE .EQ. 4) GO TO 30 00015900 C 00016000 C RESAMPLE THE TRACES (TIME TO PSEUDODEPTH CONVERSION) 00016100 C 00016200 DO 20 ID = 1, ND 00016300 VAL = TTOD(ID) / DT 00016400 IT = VAL + 1.0 00016500 IF (IT .EQ. NT) Y(ID) = Z(NT) 00016600 IF (IT .GE. NT) GO TO 30 00016700 R = VAL - (IT - 1) 00016800 C 00016900 C USING THE COEFFICIENTS COMPUTED EARLIER, FIND THE VALUE Z(IT+R) 00017000 C 00017100 CALL MWRINT (Z, NT, IT, R, LR, COEF, ANS) 00017200 Y(ID) = ANS 00017300 C 00017400 20 CONTINUE 00017500 C 00017600 30 IF (ND .EQ. NTP) GO TO 50 00017700 C 00017800 C ZERO FILL FROM END OF DATA TO NUMBER OF POINTS IN FFT 00017900 C 00018000 DO 40 IT = ND1, NTP 00018100 Y(IT) = 0.0 00018200 40 CONTINUE 00018300 C 00018400 50 CONTINUE 00018500 C 00018600 C COMPUTE FFT 00018700 C 00018800 CALL S2DFT2 (NEXPT, Y, * 8000 )00018900 C 00019000 C STORE THE FFT IN THE 2D COMPLEX STORAGE ARRAY WAV 00019100 C 00019200 CALL ARMVE (Y, C1, NTP+2) 00019300 C 00019400 DO 60 IW = 1, NWP 00019500 WAV(IKX,IW) = C1(IW) 00019600 C 00019700 60 CONTINUE 00019800 C 00019900 C END OF DO 60, THE 2D ARRAY WAV IS FULL - STORE ON DISC. 00020000 C IF DISC STORAGE IS NOT REQUIRED (IDISK = 1) BYPASS THIS STEP. 00020100 C 00020200 IF (IDISK .EQ. 1) GO TO 90 00020300 C 00020400 DO 80 IW = 1, NWP 00020500 C 00020600 C THE ARRAY IS STORED TO DISC AS TRACES ALONG THE SPACE DIMENSION 00020700 C 00020800 DO 70 IX = 1, NXBUF 00020900 C1(IX) = WAV(IX, IW) 00021000 70 CONTINUE 00021100 C 00021200 CALL FOWDSD (IOAD2, IOUTR, C1) 00021300 80 CONTINUE 00021400 C 00021500 C IF ALL OF THE TRACES HAVE BEEN PROCESSED AND STORED TO DISC, GO 00021600 C TO THE NEXT STEP. OTHERWISE, UPDATE THE INDICES AND GO BACK 00021700 C AND EXECUTE 'DO 60' AGAIN. 00021800 C 00021900 IF (NTX .EQ. NX) GO TO 90 00022000 IKX = 0 00022100 KX = KX + NXBUF 00022200 NTX = NTX + NXBUF 00022300 IF (NX .GT. NTX) GO TO 10 00022400 NTX = NX 00022500 GO TO 10 00022600 C 00022700 C ALL FFTS IN TIME HAVE BEEN COMPLETED. THE NEXT STEP IS TO COMPUTE 00022800 C THE FFTS IN SPACE. THE TRACES ALONG THE SPACE DIMENSION HAVE BEEN 00022900 C BROKEN INTO PIECES OF LENGTH NXBUF. THERE ARE NWP TRACES IN THE 00023000 C SPACE DIMENSION WHICH MUST BE TRANSFORMED. 00023100 C 00023200 C NGRI IS THE NUMBER OF NXBUF LENGTHS REQUIRED FOR NX SAMPLES. 00023300 C NGRO IS THE NUMBER OF NXBUF LENGTHS REQUIRED FOR NXP SAMPLES. 00023400 C 00023500 90 NGRI = 1 + ((NX - 1) / NXBUF) 00023600 NGRO = 1 + ((NXP - 1) / NXBUF) 00023700 C 00023800 C FOURIER TRANSFORM, X--->KX FOR ALL NWP TRACES IN THE SPACE DIMEN.00023900 C FOR EACH TRACE IN THE SPACE DIMENSION, THERE MUST BE NX/NXBUF 00024000 C RECORDS READ FROM DISC AND NXP/NXBUF RECORDS WRITTEN TO DISC. 00024100 C 00024200 DO 190 IINR = 1, NWP 00024300 C 00024400 C IF THE 2D FFT FITS IN THE 2D MEMORY BUFFER, THE DATA IS ALREADY IN00024500 C MEMORY AND DOES NOT NEED TO BE READ FROM DISC. IF THIS IS THE 00024600 C CASE MOVE THE TRACE IN THE SPACE DIMENSION FROM WAV TO C1. 00024700 C OTHERWISE READ THE TRACE FROM DISC. 00024800 C 00024900 IF (IDISK .NE. 1) GO TO 110 00025000 C 00025100 DO 100 IX = 1, NX 00025200 C1(IX) = WAV(IX, IINR) 00025300 100 CONTINUE 00025400 C 00025500 GO TO 130 00025600 C 00025700 C READ THE TRACE IN SPACE FROM THE DISC. 00025800 C STORE THE COMPLEX TRACE IN THE C1 ARRAY. 00025900 C 00026000 110 INCP = IINR 00026100 C 00026200 C NGRI SEGMENTS MUST BE READ FROM DISC TO OBTAIN THE ENTIRE TRACE. 00026300 C 00026400 DO 120 IG = 1, NGRI 00026500 INDX = 1 + (IG - 1) * NXBUF 00026600 CALL FORDSD (IOAD2, INCP, C1(INDX)) 00026700 INCP = IINR + IG * NWP 00026800 C 00026900 120 CONTINUE 00027000 C 00027100 C ZERO FILL THE TRACE OUT TO THE NUMBER OF POINTS IN THE FFT. 00027200 C 00027300 130 IF (NX1 .GT. NXP) GO TO 150 00027400 C 00027500 DO 140 IX = NX1, NXP 00027600 140 C1(IX) = CMPLX(0.0, 0.0) 00027700 C 00027800 C COMPUTE THE FFT IN SPACE. 00027900 C 00028000 150 CALL S2FFT2 (NEXPX, C1, * 8000 )00028100 C 00028200 C IF THE ENTIRE 2D FFT FITS IN THE MEMORY BUFFER WAV, RETURN THE 00028300 C TRANSFORMED SPACE TRACE TO WAV. OTHERWISE WRITE THE TRANSFORMED 00028400 C SPACE TRACE BACK TO DISC IN THE SAME LOCATIONS FROM WHICH IT WAS 00028500 C READ. 00028600 C 00028700 IF (IDISK .EQ. 1) GO TO 170 00028800 C 00028900 C WRITE THE TRANSFORMED SPACE TRACE TO DISC. NGRO SEGMENTS OF 00029000 C LENGTH NXBUF ARE REQUIRED TO BE PUT ON DISC. 00029100 C 00029200 INCP = IINR 00029300 C 00029400 DO 160 IG = 1, NGRO 00029500 INDX = 1 + (IG - 1) * NXBUF 00029600 CALL FOWDSD (IOAD2, INCP, C1(INDX)) 00029700 INCP = IINR + IG * NWP 00029800 C 00029900 160 CONTINUE 00030000 C 00030100 GO TO 190 00030200 C 00030300 C STORE THE TRANSFORMED SPACE TRACE IN THE WAV ARRAY. 00030400 C 00030500 170 DO 180 IX = 1, NXP 00030600 WAV(IX, IINR) = C1(IX) 00030700 180 CONTINUE 00030800 C 00030900 190 CONTINUE 00031000 C F-K MAPPING STEP *********************************************** 00031100 C 00031200 C BECAUSE ONLY DIPS OF LESS THAN OR EQUAL TO 90 DEGREES PROPAGATE, 00031300 C THE 90 DEGREE DIP LINE MUST BE DEFINED. THE POINT ALONG THE 00031400 C SPATIAL FREQUENCY AXIS WHERE THIS LINE CROSSES THE MAXIMUM TEMP- 00031500 C ORAL FREQUENCY IS COMPUTED AND ALL SPATIAL FREQUENCIES ABOVE THIS 00031600 C PONT ARE ZEROED OUT. THE PARAMETER MAXP IS THE INDEX FOR 00031700 C THE FIRST POSITIVE KX VALUE TO BE ZEROED. THE PARAMETER MAXN 00031800 C IS THE INDEX FOR THE FIRST NEGATIVE KX VALUE TO BE ZEROED. 00031900 C IN MOST DATA, NO SPATIAL FREQUENCIES WILL BE ZEROED 00032000 C BECAUSE THE 90 DEGREE DIP LINE WILL CROSS THE MAXIMUM TEMPORAL 00032100 C FREQUENCY AT A KX VALUE GREATER THAN THE NYQUIST IN KX (NYQX). 00032200 C THIS MEANS THAT IN MOST DATA, SOME SPATIAL ALIASING EXISTS. 00032300 C HOWEVER, MOST DATA DOES NOT CONTAIN DIPS TO 90 DEGREES SO THERE 00032400 C WILL BE LITTLE ALIASED DATA PRESENT. 00032500 C 00032600 NYQX = NXP / 2 + 1 00032700 MAX = INT(1.0 + (NWP - 1.0) * DKZ / DKX) 00032802 MAXP = MAX + 1 00032900 MAXN = NXP - MAXP + 2 00033000 IF (MAXN .LT. MAXP) GO TO 200 00033100 WRITE (KPPRNT, 9000) NYQX, MAXP, MAXN 00033200 200 CONTINUE 00033300 C 00033400 C GENERATE INTERPOLATION COEFFICIENTS FOR F-K MAPPING 00033500 C 00033600 LR = 12 00033700 BW = 0.85 00033800 CALL MCOFGN (LR, BW, RW, G, FF, A, COEF) 00033900 C 00034000 C THE F-K MAPPING OCCURS ONLY ALONG TRACES OF CONSTANT KX. THAT IS 00034100 C ONE TEMPORAL FREQUENCY AT A GIVEN KX IS MAPPED TO ANOTHER TEMPORAL00034200 C FREQUENCY AT THE SAME KX. NO MOVEMENT OF ENERGY OCCURS ALONG THE 00034300 C KX DIMENSION. IN ORDER TO PERFORM THIS MAPPING STEP IT IS 00034400 C NECESSARY TO HAVE TRACES OF CONSTANT KX SCANNING THE TEMPORAL 00034500 C FREQUENCY DIMENSION. THIS IS CARRIED OUT BY READING THE DATA 00034600 C FROM DISC IN GROUPS OF MEMORY BUFFERS. THE DATA FOR ALL TEMPORAL 00034700 C FREQUENCIES AND FOR NXBUF SPATIAL FREQUENCIES ARE READ IN AND 00034800 C MAPPED. AFTER MAPPING, THIS MEMORY BUFFER IS WRITTEN TO DISC 00034900 C IN THE SAME LOCATIONS FROM WHICH IT WAS READ AND THE NEXT MEMORY 00035000 C BUFFER IS READ IN AND MAPPED. 00035100 C 00035200 INCP = 1 00035300 INCP2 = 1 00035400 C 00035500 C DO 330 FOR EACH MEMORY BUFFER. IF NO DATA HAS BEEN STORED ON 00035600 C DISC NGRO WILL BE EQUAL TO 1. 00035700 C 00035800 DO 330 IG = 1, NGRO 00035900 C 00036000 C IF THE ENTIRE SECTION IS IN MEMORY DO NOT READ FROM DISC. 00036100 C 00036200 IF (IDISK .EQ. 1) GO TO 230 00036300 C 00036400 DO 220 IW = 1, NWP 00036500 CALL FORDSD (IOAD2, INCP, C1) 00036600 C 00036700 DO 210 IX = 1, NXBUF 00036800 210 WAV(IX, IW) = C1(IX) 00036900 C 00037000 220 CONTINUE 00037100 C 00037200 230 CONTINUE 00037300 C 00037400 C COMPUTE THE PARAMETERS FOR INDEXING THE WAV ARRAY AND FOR 00037500 C COMPUTING THE ABSOLUTE VALUE OF THE SPATIAL FREQUENCY IN 00037600 C CYCLES / FOOT. IXB AND IXE ARE POINTERS WHICH EVENTUALLY SCAN 00037700 C FROM 1 TO NXP. IKX IS A POINTER WHICH SCANS FROM 1 TO NXBUF FOR 00037800 C EACH MEMORY BUFFER TO BE PROCESSED. 00037900 C NO MAPPING OCCURS FOR ZERO CYCLES/FOOT SO IXB STARTS AT 2 INSTEAD 00038000 C OF 1. 00038100 C 00038200 IXB = 1 + (IG - 1) * NXBUF 00038300 IF (IXB .EQ. 1) IXB = 2 00038400 IXE = IG * NXBUF 00038500 IF (IXE .GT. NXP) IXE = NXP 00038600 IKX = 1 00038700 IF (IXB .EQ. 2) IKX = 2 00038800 C 00038900 C FOR EACH MEMORY BUFFER 'DO 300' 00039000 C 00039100 DO 300 IX = IXB, IXE 00039200 C 00039300 C IMX IS THE INDEX USED TO COMPUTE THE ABSOLUTE SPATIAL FREQUENCY. 00039400 C WHILE IX GOES FROM 2 TO NXP, IMX GOES FROM 2 TO NYQX BACK TO 2. 00039500 C 00039600 IMX = IX 00039700 IF (IX .GT. NYQX) IMX = 2 * NYQX - IX 00039800 C 00039900 C IF IX IS GREATER THAN MAXP AND LESS THAN MAXN THE SPATIAL 00040000 C FREQUENCY LIES IN THE ZONE OF NO PROPAGATION AND IS ZEROED OUT 00040100 C (GO TO 270). 00040200 C 00040300 IF (IX .GE. MAXP .AND. IX .LE. MAXN) GO TO 270 00040400 C 00040500 C IN ORDER TO SIMPLIFY THE INTERPOLATION IN THE F-K DOMAIN 00040600 C THE COMPLEX DATA ARE SORTED INTO REAL AND IMAGINARY PARTS WHICH 00040700 C ARE INTERPOLATED INDEPENDENTLY. 00040800 C 00040900 DO 240 IKZ = 1, NWP 00041000 Y(IKZ) = REAL (WAV(IKX, IKZ)) 00041100 240 Z(IKZ) = AIMAG(WAV(IKX, IKZ)) 00041200 C 00041300 C AKX IS THE SPATIAL FREQUENCY IN CYCLES/FOOT. 00041400 C 00041500 AKX = (IMX - 1.0) * DKX 00041600 AKX2 = AKX * AKX 00041700 THFR = AKX * COS(DIPS) / SIN(DIPS) 00041800 THMN = AKX * COS(DIPW) / SIN(DIPW) 00041900 AKZDIF = THFR - THMN 00042000 C 00042100 DO 260 IKZ = 1, NWP 00042200 C 00042300 C AKZ IS THE TEMPORAL FREQUENCY (IN CYCLES/SEC) FOR THE OUTPUT 00042400 C F-K SECTION. F IS THE TEMPORAL FREQUENCY IN CYCLES/SEC OF 00042500 C THE INPUT F-K SECTION. THE VALUE STORED AT F IS TO BE MOVED 00042600 C TO THE LOCATION AKZ. SINCE THE VALUE OF F DOES NOT LIE ON AN 00042700 C EXISTING SAMPLE IT MUST BE INTERPOLATED. IW IS THE INDEX TO THE 00042800 C TEMPORAL FREQUENCY AXIS SUCH THAT F LIES BETWEEN IW AND IW+1. 00042900 C R IS THE FRACTIONAL DISTANCT BETWEEN IW AND IW+1. 00043000 C SCALF IS THE DIP FILTER SCALING FACTOR. 00043100 C 00043200 AKZ = (IKZ - 1.0) * DKZ 00043300 SCALF = 1.0 00043400 DIF = THFR - AKZ 00043500 IF (DIF .LT. 0.0) GO TO 248 00043600 IF (AKZDIF .NE. 0.0) GO TO 245 00043700 SCALF = 0.0 00043800 GO TO 248 00043900 C 00044000 245 SCALF = SCALF - (DIF / AKZDIF) 00044100 IF (SCALF .LT. 0.0) SCALF = 0.0 00044200 C 00044300 248 F = VC * SQRT(AKX2 + AKZ * AKZ) / 2.0 00044400 F1 = F / DKD 00044500 IW = F1 + 1.0 00044600 R = F1 - (IW - 1.0) 00044700 C 00044800 C IF IW IS GREATER THAN THE MAXIMUM TEMPORAL FREQUENCY AVAILABLE, 00044900 C THE DESIRED VALUE IS NOT AVAILABLE SO IT IS ZEROED OUT. 00045000 C 00045100 IF (IW .LT. NWP) GO TO 250 00045200 WAV(IKX, IKZ) = CMPLX(0.0, 0.0) 00045300 GO TO 260 00045400 C 00045500 C WIENER INTERPOLATION TO OBTAIN THE VALUE TO STORE AT AKZ. 00045600 C 00045700 250 CALL MWRINT (Y, NWP, IW, R, LR, COEF, WR) 00045800 CALL MWRINT (Z, NWP, IW, R, LR, COEF, WI) 00045900 WA = CMPLX(WR, WI) 00046000 C 00046100 C OBLIQ IS THE AMPLITUDE SCALE FACTOR IN THE MAPPING. 00046200 C 00046300 OBLIQ = VC * AKZ / (F * 2.0) 00046400 WAV(IKX, IKZ) = WA * OBLIQ * SCALF 00046500 260 CONTINUE 00046600 C 00046700 GO TO 290 00046800 C 00046900 C IF THE SPATIAL FREQUENCY IS GREATER THAN MAXP AND LESS THAN MAXN, 00047000 C ZERO OUT ALL OF THE TEMPORAL FREQUENCIES. 00047100 C 00047200 270 DO 280 I = 1, NWP 00047300 C 00047400 280 WAV(IKX, I) = CMPLX(0.0, 0.0) 00047500 C 00047600 290 IKX = IKX + 1 00047700 C 00047800 300 CONTINUE 00047900 C 00048000 C AFTER PERFORMING THE MAPPING FOR ALL SPATIAL FREQUENCIES WRITE 00048100 C THE MEMORY BUFFER TO DISC. 00048200 C 00048300 C IF THE ENTIRE 2D FFT FITS IN MEMORY DO NOT WRITE TO DISC. 00048400 C 00048500 IF (IDISK .EQ. 1) GO TO 340 00048600 C 00048700 DO 320 IW = 1, NWP 00048800 C 00048900 DO 310 IX = 1, NXBUF 00049000 310 C1(IX) = WAV(IX, IW) 00049100 C 00049200 CALL FOWDSD (IOAD2, INCP2, C1) 00049300 320 CONTINUE 00049400 C 00049500 330 CONTINUE 00049600 C 00049700 340 CONTINUE 00049800 C 00049900 C INVERSE 2D FFT ************************************************* 00050000 C 00050100 C THE MAPPING STEP IS COMPLETE. NOW THE INVERSE 2D FFT MUST BE 00050200 C CARRIED OUT. FIRST DO THE INVERSE FFTS IN SPACE AND THEN IN 00050300 C IN TIME. THERE ARE NWP SPACE FFTS TO BE PERFORMED. 00050400 C 00050500 C INVERSE FOURIER TRANSFORM, KX --> X 00050600 C 00050700 DO 420 IINR = 1, NWP 00050800 C 00050900 C IF THE ENTIRE 2D FFT IS IN MEMORY DO NOT READ FROM DISC. STORE 00051000 C THE SPATIAL FREQUENCIES FOR EACH TEMPORAL FREQUENCY IN C1. 00051100 C 00051200 IF (IDISK .NE. 1) GO TO 360 00051300 C 00051400 DO 350 IX = 1, NXP 00051500 C1(IX) = WAV(IX, IINR) 00051600 C 00051700 350 CONTINUE 00051800 C 00051900 GO TO 380 00052000 C 00052100 C NGRO SEGMENTS MUST BE READ FROM DISC FOR EACH TRACE IN THE SPATIAL00052200 C FREQUENCY DIMENSION. 00052300 C 00052400 360 INCP = IINR 00052500 C 00052600 DO 370 IG = 1, NGRO 00052700 INDX = 1 + (IG - 1) * NXBUF 00052800 CALL FORDSD (IOAD2, INCP, C1(INDX)) 00052900 INCP = IINR + IG * NWP 00053000 C 00053100 370 CONTINUE 00053200 C 00053300 380 CONTINUE 00053400 C 00053500 C HAVING STORED THE SPATIAL FREQENCY TRACE IN C1. DO THE FFT. 00053600 C 00053700 CALL S2FFI2 (NEXPX, C1, * 8000 )00053800 C 00053900 C IF THE ENTIRE 2D FFT IS IN MEMORY STORE THE RESULT BACK IN WAV. 00054000 C 00054100 IF (IDISK .NE. 1) GO TO 400 00054200 C 00054300 DO 390 IX = 1, NXP 00054400 390 WAV(IX, IINR) = C1(IX) 00054500 C 00054600 GO TO 420 00054700 C 00054800 C WRITE THE RESULT OF THE FFT BACK TO DISC IN THE SAME LOCATIONS 00054900 C FROM WHICH THE INPUT WAS READ. BECAUSE ONLY NX VALUES NEED BE 00055000 C RETAINED, NGRI SEGMENTS ARE REQUIRED TO BE WRITTEN. 00055100 C 00055200 400 INCP = IINR 00055300 C 00055400 DO 410 IG = 1, NGRI 00055500 INDX = 1 + (IG - 1) * NXBUF 00055600 CALL FOWDSD (IOAD2, INCP, C1(INDX)) 00055700 INCP = IINR + IG * NWP 00055800 C 00055900 410 CONTINUE 00056000 C 00056100 420 CONTINUE 00056200 C 00056300 C THE SPATIAL INVERSE FFTS ARE COMPLETE. NOW DO THE TEMPORAL FFTS. 00056400 C BECAUSE THE OUTPUT OF THE FFT IS THE DESIRED TRACE, COMPLETE 00056500 C THE MIGRATION PROCESS BY MAPPING THE TRACE BACK TO THE ORIGINAL 00056600 C SAMPLING AND STORE THE RESULT ON THE TRACE SAVE FILE. 00056700 C THE INTERPOLATION COEFFICIENTS FOR MAPPING THE TRACES ARE THE SAME00056800 C ONES COMPUTED FOR THE F-K MAPPING STEP 00056900 C IN ORDER TO PERFORM THE FFTS IN TIME, THE DATA ARE READ FROM DISC 00057000 C ONE MEMORY BUFFER AT A TIME. THERE ARE NGRI MEMORY BUFFERS TO BE 00057100 C PROCESSED. THE FFTS ARE PERFORMED ON EACH MEMORY BUFFER BEFORE 00057200 C THE NEXT BUFFER IS READ IN. 00057300 C WE HAVE RETAINED ONLY THE POSITIVE TEMPORAL FREQUENCIES THROUGH 00057400 C THE MIGRATION PROCESS. HOWEVER, BECAUSE THE INPUT DATA WERE REAL 00057500 C (NOT COMPLEX), THE REAL PART OF THE TEMPORAL FREQUENCIES IS 00057600 C SYMETRICAL ABOUT ZERO HZ AND THE IMAGINARY PART IS ANTISYMETRICAL.00057700 C THEREFORE THE NEGATIVE FREQUENCIES NEEDED FOR THE INVERSE FFT IN 00057800 C TIME CAN BE GENERATED. 00057900 C 00058000 INCP = 1 00058100 IKX = 1 00058200 NWPM = NWP - 1 00058300 C 00058400 C READ ALL OF THE TEMPORAL FREQUENCIES FROM DISC FOR NXBUF TRACES. 00058500 C NGRI GROUPS OF NXBUF TRACES MUST BE READ FROM DISC. 00058600 C 00058700 DO 490 IG = 1, NGRI 00058800 C 00058900 C IF THE ENTIRE 2D FFT IS IN MEMORY DO NOT READ FROM DISC. 00059000 C 00059100 IF (IDISK .EQ. 1) GO TO 450 00059200 C 00059300 DO 440 IW = 1, NWP 00059400 CALL FORDSD (IOAD2, INCP, C1) 00059500 C 00059600 DO 430 IX = 1, NXBUF 00059700 430 WAV(IX, IW) = C1(IX) 00059800 C 00059900 440 CONTINUE 00060000 C 00060100 C DO THE INVERSE FFT AND MAPPING FOR EACH OF THE TRACES IN 00060200 C THE MEMORY BUFFER. 00060300 C 00060400 450 DO 480 IX = 1, NXBUF 00060500 C 00060600 C MOVE ALL OF THE POSITIVE FREQUENCIES FOR A GIVEN TRACE TO THE 00060700 C C1 ARRAY. FILL IN THE NEGATIVE FREQUENCIES BY TAKING THE 00060800 C COMPLEX CONJUGATE OF THE POSITIVE FREQUENCIES. 00060900 C 00061000 DO 460 IW = 2, NWPM 00061100 C1(IW) = WAV(IX, IW) 00061200 JW = NTP - IW + 2 00061300 C1(JW) = CONJG (C1(IW)) 00061400 C 00061500 460 CONTINUE 00061600 C 00061700 VAL = REAL (WAV(IX, 1)) 00061800 C1(1) = CMPLX (VAL, 0.0) 00061900 VAL = REAL (WAV(IX, NWP)) 00062000 C1(NWP) = CMPLX (VAL, 0.0) 00062100 C 00062200 C INVERSE FFT IN TIME. 00062300 C 00062400 CALL ARMVE (C1, Y, NTP+2) 00062500 CALL S2DFI2 (NEXPT, Y, * 8000 )00062600 C 00062700 C MAP THE OUTPUT TRACE FROM ND POINTS TO NT POINTS ACCORDING TO 00062800 C THE DATA IN ARRAY TTOD. 00062900 C 00063000 IF (MODE .EQ. 3) GO TO 475 00063100 C 00063200 DO 470 IT = 1, NT 00063300 VAL = DTOT(IT) / DD 00063400 ID = INT(VAL) + 1 00063500 R = VAL - (ID - 1.0) 00063600 C 00063700 C WIENER INTERPOLATION TO OBTAIN THE OUTPUT VALUE. 00063800 C 00063900 CALL MWRINT (Y, ND, ID, R, LR, COEF, ANS) 00064000 Z(IT) = ANS 00064100 470 CONTINUE 00064200 C 00064300 C OUTPUT THE MIGRATED TRACE 00064400 C 00064500 IF (MODE .EQ. 1) IOREC = ((ILINE-1) * CDPFLD + JF - 1) * NXM + IKX00064600 IF (MODE .EQ. 2) IOREC = ((IKX-1) * CDPFLD + JF - 1) * NXM + ILINE00064700 IF (MODE .EQ. 4) IOREC = ((IKX-1) * CDPFLD + JF - 1) * NXM + ILINE00064800 475 IF (MODE .EQ. 3) IOREC = ((ILINE-1) * CDPFLD + JF - 1) * NXM + IKX00064900 CALL FORDSD (IOADR, IOREC, WBUF) 00065000 C 00065100 CALL ARSET (WBUF(THL+1), NT, 0.0) 00065200 IF (MODE .NE. 3) CALL ARMVE (Z, WBUF(THL+1), NT) 00065300 IF (MODE .EQ. 3) CALL ARMVE (Y, WBUF(THL+1), ND) 00065400 C 00065500 IOREC = IOREC - 1 00065600 C 00065700 CALL FOWDSD (IOADR, IOREC, WBUF) 00065800 C 00065900 IKX = IKX + 1 00066000 IF (IKX .GT. NX) GO TO 500 00066100 C 00066200 480 CONTINUE 00066300 C 00066400 490 CONTINUE 00066500 C 00066600 500 RETURN 00066700 C 00066800 8000 IER = -1 00066900 GO TO 500 00067000 C 00067100 9000 FORMAT (3X,'NYQUIST INDEX = ',I10,', ZERO OUT ABOVE ',I10, 00067200 * ', AND BELOW ',I10) 00067300 END 00067400