CTITLESAHCI _ COMPUTE AVO FOR PROCESS AVEL C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** C@PROCESS VECTOR(LEVEL(2) REPORT(XLIST)) C CA AUTHOR H. W. SWAN CA DESIGNER H. W. SWAN CA SYSTEM IBM/CRAY CA LANGUAGE VS FORTRAN VERSION 2.2 CA WRITTEN 12-03-87 CA LAST REVISED 08-17-88 C REVISED 04-26-88 HWS CHANGED SCOND TO MACONV C 06-13-88 HWS ADDED PHASE ROTATION C 08-17-88 HWS ALLOWS NEGATIVE V0 C 11-21-89 HWS PARM(14)= A-B ROTATION C ANGLE (10'THS OF DEGREES) C REVISED 05-23-91 JJC CHANGED IMPLICIT NONE TO (A-Z). C REVISED 12-17-91 JJC MODIFIED TO MEET SPARC STANDARDS. CA CA CALLING SEQUENCE: CA CALL SAHCI(S, H1, H2, SLOTH, SLOTHI, PARM, AOUT, BOUT, CA MAXSMP, MAXH, NPARM, NSAMP, NH, T0, TSAMP, AVOPCT) CA CA IN/OUT ARGUMENT TYPE DESCRIPTION CA ------ -------- ---- ----------- CA CA IN/OUT S (R4) COLLECTED WEIGHTED SUMS CA (2-D ARRAY, DIMENSIONED MAXSMP X 8) CA CA INPUT: S(1,1) = S(1,SY) = SUM { D(X) } CA S(1,2) = S(1,SX2Y) = SUM { D(X) * X**2 } CA S(1,3) = S(1,SX4Y) = SUM { D(X) * X**4 } CA S(1,4) = S(1,SX0) = SUM { X**0 } CA S(1,5) = S(1,SX2) = SUM { X**2 } CA S(1,6) = S(1,SX4) = SUM { X**4 } CA S(1,7) = S(1,SX6) = SUM { X**6 } CA S(1,8) = S(1,SX8) = SUM { X**8 } CA CA OUTPUT: ZEROED, OR USED AS SCRATCH BY THE ROUTINE CA (2-D ARRAY, DIMENSIONED MAXSMP X 8) CA CA IN H1 (R4) THE FIRST SET OF INVERSION FILTER COEFFICIENTS CA (1-D ARRAY, DIMENSIONED -MAXH:MAXH) CA CA IN H2 (R4) THE SECOND SET OF INVERSION FILTER COEFFICIENTS CA (1-D ARRAY, DIMENSIONED -MAXH:MAXH) CA CA IN SLOTH (R4) THE STACKING SLOTH WHICH WAS USED. CA `SLOTH' IS DEFINED TO BE { VELOCITY ** (-2)}. CA (1-D ARRAY, DIMENSIONED 0:NSAMP+1) CA CA IN SLOTHI (R4) THE CORRESPONDING INTERVAL SLOTH. IF NOT CA AVAILABLE, THE STACKING SLOTH MAY BE USED. CA (1-D ARRAY, DIMENSIONED 0:NSAMP+1) CA CA IN PARM (I4) AN INTEGER PARAMETER ARRAY: CA CA PARM(1) = THE HILBERT TRANSFORM FILTER LENGT CA CA (2) = THE "A" ORDER: CA 0 = UNITY CA 1 = CONVENTIONAL STACK CA 2 = ZERO-OFFSET STACK, QUADRATIC FIT CA 3 = ZERO-OFFSET STACK, QUARTIC FIT CA CA PARM(3) = THE "B" ORDER: CA 0 = UNITY CA 1 = QUADRATIC SLOPE STACK, QUADRATIC FIT CA 2 = QUADRATIC SLOPE STACK, QUARTIC FIT CA 3 = QUARTIC SLOPE CA CA PARM(4) = FLAG FOR NMO STRETCH CORRECTION: CA 0 = DO NOT CORRECT FOR NMO STRETCH CA 1 = DO CORRECT FOR NMO STRETCH CA CA PARM(5) = HORIZONTAL SPHERICAL DIVERGENCE FL CA 0 = DO NOT CORRECT FOR HORIZONTAL CA SPHERICAL DIVERGENCE EFFECTS CA 1 = DO CORRECT FOR HORIZONTAL CA SPHERICAL DIVERGENCE EFFECTS CA CA PARM(6) = VERTICAL SPHERICAL DIVERGENCE FLAG CA 0 = DO NOT CORRECT FOR VERTICAL SPHERICA CA DIVERGENCE EFFECTS CA 1 = DO CORRECT FOR VERTICAL SPHERICAL CA DIVERGENCE EFFECTS CA CA PARM(7) = V0 FOR SPHERICAL DIVERGENCE CA CORRECTION. CA CA PARM(8) = TYPE OF OUTPUT: CA 0: AOUT = RE{AB*}, BOUT = IM{AB*}; CA 1: AOUT = RE{A}, BOUT = IM{AB*}; CA 2: AOUT = RE{A}, BOUT = RE{AB*}; CA 3: AOUT = RE{A}, BOUT = RE{B}. CA CA PARM(9) = THE LENGTH OF THE DATA FILTER CA CA PARM(10): AGC WINDOW LENGTH FOR |A| CA SCALING. (IF 0, THEN NO AGC.) CA (IF <0, THEN SCALE BY |A|**2) CA CA PARM(11): AGC WINDOW LENGTH FOR |B| CA SCALING. (IF 0, THEN NO AGC.) CA (IF <0, THEN SCALE BY |B|**2) CA CA PARM(12): TAPER LENGTH IN SAMPLES AT CA BOTH EDGES OF THE ANALYSIS WINDOW CA CA PARM(13): THE PHASE ANGLE BY WHICH TO CA ROTATE THE COMPUTES SLOPES AND CA INTERCEPTS. CA CA PARM(14): THE ROTATION ANGLE IN THE A-B CA PLANE, IN TENTHS OF DEGREES. CA CA OUT AOUT (R4) THE FIRST HYDROCARBON INDICATOR ARRAY CA (1-D ARRAY, DIMENSIONED MAXSMP) CA CA OUT BOUT (R4) THE SECOND HYDROCARBON INDICATOR ARRAY CA (1-D ARRAY, DIMENSIONED MAXSMP) CA CA IN MAXSMP I4 THE DIMENSIONED SIZES OF THE S, SLOTH, SLOTHI, CA AOUT, AND BOUT ARRAYS CA CA IN MAXH I4 THE SIZE OF THE H ARRAY CA CA IN NPARM I4 THE SIZE OF THE PARM ARRAY CA CA IN NSAMP I4 THE ACTUAL SIZES OF THE S, SLOTH, SLOTHI, AOUT, CA AND BOUT ARRAYS CA CA IN NH I4 THE ACTUAL SIZE OF THE H ARRAY CA CA IN T0 R4 THE TIME IN SECONDS OF THE FIRST GIVEN SAMPLE CA CA IN TSAMP R4 THE SAMPLING INTERVAL IN SECONDS CA CA IN AVOPCT R4 THE AVO STABILIZATION TERM, USED IN CA MATRIX INVERSION. CA CA CA PURPOSE: CA CA THIS SUBROUTINE COMPUTES VARIOUS HYDROCARBON INDICATORS FROM CA COLLECTED WEIGHTED SUMS. CA CA C C*********************************************************************** C*** **** C*** SUBROUTINES CALLED BY 'SAHCI': **** C*** **** C*********************************************************************** C C C ROUTINES CONTAINED AS PART OF THE 'AVEL' PACKAGE: C ------------------------------------------------- C C MACONV - PERFORMS A DIRECT CONVOLUTION C MAHLBT - PERFORMS A HILBERT TRANSFORM OF THE DATA C MASMTH - SMOOTHS THE DATA WITH A RECURSIVE BOXCAR FILTER C C ESSL/SCILIB ROUTINES USED: C -------------------------- C C SCOPY - COPIES ONE VECTOR INTO ANOTHER C C C*********************************************************************** C*** **** C*** CALLING PARAMETERS OF 'SAHCI': **** C*** **** C*********************************************************************** C C SUBROUTINE SAHCI(S, H1, H2, SLOTH, SLOTHI, PARM, AOUT, BOUT, * MAXSMP, MAXH, NPARM, NSAMP, NH, T0, TSAMP, AVOPCT) C CJJ IMPLICIT NONE IMPLICIT INTEGER (A-Z) C EXTERNAL SCOPY, MACONV, MAHLBT, MASMTH INTEGER SY, SX2Y, SX4Y, SX0, SX2, SX4, SX6, SX8 PARAMETER (SY=1, SX2Y=2, SX4Y=3, SX0=4, SX2=5, SX4=6) PARAMETER (SX6=7, SX8=8) INTEGER MAXSMP, MAXH, NPARM, NH, AINDEX, BINDEX, PHROT INTEGER PARM(NPARM), DSTRCH, SPHFLG, FTRLEN, VFLAG INTEGER TNDX(6,3), DNDX(6), ORDER(2), KEY(2, 0:3) INTEGER IHILB, IDSMTH, ISTG, JA, J, NSAMP, MNDX INTEGER N0, N2, N4, ND, I, IAFLG, IBFLG, TAPER DOUBLE PRECISION P, Q, DSDT, DF, EP, AVO, AVOP DOUBLE PRECISION S0, S2, S4, S6, S8, S4P, AN REAL S, H1, H2, SLOTH, SLOTHI, AOUT, BOUT, T0, TSAMP REAL AVOPCT, V0, PFACT, DFACT, TIME, PI, TEMP REAL SINPR, COSPR, SV2, ROTAB DIMENSION S(MAXSMP,8), H1(-MAXH:MAXH), H2(-MAXH:MAXH), * SLOTH(0:MAXSMP+1), SLOTHI(0:MAXSMP+1), AOUT(MAXSMP), * BOUT(MAXSMP) DOUBLE PRECISION T(10), D(3), DENOM, EPS PARAMETER (EPS = 1D-40, PI=3.141593) LOGICAL OKCORR DATA TNDX/ 2, 4, 10, 3, 9, 7, * 1, 3, 9, 2, 8, 6, * 1, 1, 7, 1, 6, 5 / DATA DNDX/ 1, 2, 3, 2, 3, 3 / DATA KEY / 6, 7, 1, 7, 1, 6, 1, 2 / C C FIGURE OUT WHAT TO COMPUTE. C FTRLEN = PARM(1) ORDER(1) = PARM(2) ORDER(2) = PARM(3) DSTRCH = PARM(4) SPHFLG = PARM(5) VFLAG = PARM(6) V0 = PARM(7) IHILB = PARM(8) IDSMTH = PARM(9) IAFLG = PARM(10) IBFLG = PARM(11) TAPER = PARM(12) PHROT = PARM(13) ROTAB = PARM(14) / 10.0 OKCORR = ORDER(1) .NE. 0 .AND. ORDER(2) .NE. 0 C C ADD WHITE NOISE TO SUM(X(I)**N), FOR N = 0, 4, 8. C AVO = 10.0 ** (-AVOPCT / 20.0) AVOP= AVO / (1D0 + AVO) DO 50 ISTG = 8, 4, -2 DO 50 I = 1,NSAMP AN = S(I, SX0) * AVO 50 S(I, ISTG) = S(I, ISTG) + AN C C COMPUTE THE RAW ANSWERS IN 'AOUT' AND 'BOUT'. C DO 200 ISTG = 1,2 JA = ORDER(ISTG) IF(JA .LT. 0 .OR. JA .GT. 3) JA = 0 IF(JA .EQ. 0) THEN DO 80 J=1, NSAMP 80 BOUT(J) = 1.0 GO TO 200 ENDIF MNDX = (ISTG-1)*3 + JA N0 = TNDX(MNDX, 1) N2 = TNDX(MNDX, 2) N4 = TNDX(MNDX, 3) ND = DNDX(MNDX) T(1) = 0.0 DO 100 I=1, NSAMP S0 = S(I, SX0) S2 = S(I, SX2) S4 = S(I, SX4) S6 = S(I, SX6) S8 = S(I, SX8) EP = S4 * AVOP S4P = S4 - EP C T(2) = S0 T(3) =-S2 T(4) = S4 T(5) = S0*S4 - S2*S2 T(6) = S2*S4P- S0*S6 T(7) = S2*S6 - S4*S4 + S4*EP T(8) = S0*S8 - S4P*S4P T(9) = S6*S4P- S2*S8 T(10)= S4*S8 - S6*S6 D(1) = T(2) D(2) = T(5) D(3) = S8*T(5) + S6*T(6) + S4*T(7) * + EP*(S4*S4P - S2*S6) DENOM = D(ND) BOUT(I) = 0.0 IF(DENOM .GT. EPS) * BOUT(I) = (T(N0)*S(I, SY) + T(N2)*S(I, SX2Y) + * T(N4)*S(I, SX4Y)) / DENOM 100 CONTINUE IF(ISTG .EQ. 1) CALL SCOPY(NSAMP, BOUT, 1, AOUT, 1) 200 CONTINUE C C TAPER 'AOUT' AND 'BOUT' AT BOTH EDGES. C DO 240 I=1,TAPER PFACT = FLOAT(I) / TAPER IF(I .LE. NSAMP) THEN AOUT(I) = AOUT(I) * PFACT BOUT(I) = BOUT(I) * PFACT ENDIF 240 CONTINUE DO 250 I=NSAMP-TAPER+1, NSAMP PFACT = FLOAT(NSAMP - I + 1) / TAPER IF(I .GT. 0) THEN AOUT(I) = AOUT(I) * PFACT BOUT(I) = BOUT(I) * PFACT ENDIF 250 CONTINUE C C ZERO OUT UNWANTED PARTS OF S: C DO 270 J=1, 8 DO 270 I=1,MAXSMP 270 S(I,J) = 0.0 C C COPY S(1) <-- 'A'; S(2) <-- 'B', IN CASE NO FILTERING IS REQUIRED. C CALL SCOPY(NSAMP, AOUT, 1, S(1,1), 1) CALL SCOPY(NSAMP, BOUT, 1, S(1,2), 1) C C NOW FILTER BOTH AEST(T) AND BEST(T) WITH THE H1 INVERSION FILTER. C C IF(ORDER(1) .NE. 0) C * CALL SCOND(H1(-NH), 1, AOUT, 1, S, 1, C * 2*NH+1, NSAMP, NH, NSAMP) C IF(ORDER(2) .NE. 0) C * CALL SCOND(H1(-NH), 1, BOUT, 1, S(1,2), 1, C * 2*NH+1, NSAMP, NH, NSAMP) IF(ORDER(1) .NE. 0) * CALL MACONV(AOUT, H1(-NH), S(1,1), NSAMP, 2*NH+1, NH+1) IF(ORDER(2) .NE. 0) * CALL MACONV(BOUT, H1(-NH), S(1,2), NSAMP, 2*NH+1, NH+1) C C DO THE DESTRETCHING, IF CALLED FOR: C S(3) <-- A * H2(T). C IF(DSTRCH .EQ. 1 .AND. OKCORR) THEN C CALL SCOND(H2(-NH), 1, AOUT, 1, S(1,3), 1, C * 2*NH+1, NSAMP, NH, NSAMP) CALL MACONV(AOUT, H2(-NH), S(1,3), NSAMP, 2*NH+1, NH+1) ENDIF C C PERFORM IN-PLACE CORRECTION OF 'B': C PFACT = 0.0 DFACT = 0.0 IF(OKCORR) DFACT = 0.5 IF(SPHFLG .NE. 0 .AND. OKCORR) PFACT = 0.5 DO 400 I=1, NSAMP TIME = (I-1)*TSAMP + T0 SV2 = SLOTH(I) * V0**2 C P = PFACT * (2.0 - SIGN(SV2, V0)) P = PFACT * (2.0 - SV2) IF(V0 .LT. 0) * P = PFACT * (1.0 + SV2) / SV2 Q = SLOTH(I) / SLOTHI(I) DSDT = (1D0*SLOTH(I+1) - SLOTH(I-1))/TSAMP/2.0 DF = DFACT*(1.0 - TIME*DSDT/SLOTH(I)) S(I,2) = (S(I,2) + P*S(I,1) + DF*S(I,3))/Q 400 CONTINUE C C CORRECT FOR VERTICAL SPHERICAL DIVERGENCE EFFECT, IF SO REQUESTED. C IF(VFLAG .NE. 0) THEN DO 450 I=1, NSAMP TIME = (I-1)*TSAMP + T0 P = TIME / (SLOTH(I) * V0**2) S(I,1) = S(I,1) * P 450 S(I,2) = S(I,2) * P ENDIF C C NOW PERFORM A HILBERT TRANSFORM OF THE DATA. C S(4) <-- HI { S(1) }; S(5) <-- HI { S(2) }. C DO 500 I=1,2 IF(ORDER(I) .EQ. 0) GO TO 500 CALL MAHLBT(S(1,I), S(1,I+3), NSAMP, I, NSAMP, FTRLEN) 500 CONTINUE C C ROTATE THE RESULTS IN TIME, IF CALLED FOR. C IF(PHROT .NE. 0) THEN COSPR = COS(PI*PHROT/180.0) SINPR = SQRT(1.0 - COSPR**2) DO 550 I=1, NSAMP TEMP = S(I,1)*SINPR + S(I,4)*COSPR S(I,1) = S(I,1)*COSPR - S(I,4)*SINPR S(I,4) = TEMP TEMP = S(I,2)*SINPR + S(I,5)*COSPR S(I,2) = S(I,2)*COSPR - S(I,5)*SINPR S(I,5) = TEMP 550 CONTINUE ENDIF C C ROTATE THE RESULTS IN THE A-B PLANE, IF CALLED FOR. C IF(ROTAB .NE. 0) THEN COSPR = COS(PI*ROTAB/180.0) SINPR = SQRT(1.0 - COSPR**2) DO 560 I=1, NSAMP TEMP = S(I,2)*SINPR + S(I,1)*COSPR S(I,2) = S(I,2)*COSPR - S(I,1)*SINPR S(I,1) = TEMP TEMP = S(I,5)*SINPR + S(I,4)*COSPR S(I,5) = S(I,5)*COSPR - S(I,4)*SINPR S(I,4) = TEMP 560 CONTINUE ENDIF C C PLACE REAL AND IMAGINARY PARTS OF A * CONJG(B) INTO C (S(6), S(7)), RESPECTIVELY. C DO 600 I=1, NSAMP S(I,6) = S(I,1)*S(I,2) + S(I,4)*S(I,5) S(I,7) = S(I,4)*S(I,2) - S(I,1)*S(I,5) 600 CONTINUE C C SMOOTH THE PRODUCT TERMS. C DO 650 I=1,2 CALL MASMTH(NSAMP, 1, 1, S(1,I+5), S(1,8), IDSMTH) 650 CONTINUE C C NORMALIZE BY |A|, IF IAFLG ^= 0. C IF(IAFLG .NE. 0) THEN IF(IAFLG .GT. 0) THEN DO 710 I=1, NSAMP 710 S(I,3) = SQRT(ABS(S(I,1)**2 + S(I,4)**2)) ELSE DO 715 I=1, NSAMP 715 S(I,3) = S(I,1)**2 + S(I,4)**2 ENDIF CALL MASMTH(NSAMP, 1, 1, S(1,3), S(1,8), IABS(IAFLG)) DO 720 I=1, NSAMP DENOM = S(I,3) IF(DENOM .GT. EPS) THEN S(I,6) = S(I,6) / DENOM S(I,7) = S(I,7) / DENOM ENDIF 720 CONTINUE ENDIF C C NORMALIZE BY |B|, IF IBFLG ^= 0. C IF(IBFLG .NE. 0) THEN IF(IBFLG .GT. 0) THEN DO 730 I=1, NSAMP 730 S(I,4) = SQRT(ABS(S(I,2)**2 + S(I,5)**2)) ELSE DO 740 I=1, NSAMP 740 S(I,4) = S(I,2)**2 + S(I,5)**2 ENDIF CALL MASMTH(NSAMP, 1, 1, S(1,4), S(1,8), IABS(IBFLG)) DO 750 I=1, NSAMP DENOM = S(I,4) IF(DENOM .GT. EPS) THEN S(I,6) = S(I,6) / DENOM S(I,7) = S(I,7) / DENOM ENDIF 750 CONTINUE ENDIF C C NOW COPY THE SELECTED ANSWERS TO THE OUTPUT ARRAYS. C AINDEX = KEY(1,IHILB) BINDEX = KEY(2,IHILB) CALL SCOPY(NSAMP, S(1,AINDEX), 1, AOUT, 1) CALL SCOPY(NSAMP, S(1,BINDEX), 1, BOUT, 1) C C ALL DONE! C RETURN END