CTITLESAFF3DA - BUILD TABLE FOR REJECT FILTER IN F-X DOMAIN C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CABS SAFF3DA - BUILD TABLE FOR REJECT FILTER IN F-X DOMAIN C CSUBROUTINE SAFF3DA C C SUBROUTINE SAFF3DA( H,NXT,NF,DXT,XMX,SCF, C * DF,F1,F2,F3,F4,V1,V2,V3,V4,SGN,ALB,CDB,FCT ) C C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991. C C ALL RIGHTS RESERVED. NO PART OF THIS PROGRAM MAY BE PHOTOCOPIED, C REPRODUCED, OR TRANSLATED TO ANOTHER PROGRAM LANGUAGE WITHOUT THE C PRIOR CONSENT OF ATLANTIC RICHFIELD COMPANY. C C CA DESIGNER D CORRIGAN CA AUTHOR D CORRIGAN CA LANGUAGE FORTRAN 77 CA SYSTEM IBM/CRAY CA WRITTEN 01-29-92 C REVISED 01-29-92 TWO CORRECTIONS: C 1. APPLY TAPER ONLY TO END OF C H ARRAY C 2. ACCOUNT FOR CASE WHEN C FA .GT. FB C 09-29-92 DC CHECK FOR CASE A12 = 1. C OR A34 = 1. C CA CA PURPOSE OF PROGRAM: CA CA BUILD TABLE FOR IMPLEMENTATION OF A REJECT FILTER CA IN THE F-X DOMAIN CA C*********************************************************************** C C SUBROUTINES AND FUNCTIONS CALLED: C C SAFF3DG C C*********************************************************************** C CA CA ARGUMENTS: CA CA ARGUMENT TYPE I/O MEANING CA -------- ------- --- -------- CA CA H COMPLEX O ARRAY OF FILTERS (NXT,NF) CA NXT INTEGER I NUMBER OF OFFSETS PER FREQUENCY CA NF INTEGER I NUMBER OF FREQUENCIES CA DXT REAL O OFFSET INCREMENT FOR EACH FREQ. CA XMX REAL O MAXIMUM OFFSET FOR EACH FREQ. CA SCF REAL O SCALAR MULTIPLIER FOR EACH FREQ. CA DF REAL I FREQUENCY INCREMENT CA F1 REAL I LOW PASS FREQUENCY CA F2 REAL I LOW CUT FREQUENCY CA F3 REAL I HIGH CUT FREQUENCY CA F4 REAL I HIGH PASS FREQUENCY CA V1 INTEGER I LOW PASS VELOCITY CA V2 INTEGER I LOW CUT VELOCITY CA V3 INTEGER I HIGH PASS VELOCITY CA V4 INTEGER I HIGH CUT VELOCITY CA SGN INTEGER I SIGN OF VELOCITIES CA ALB REAL I HIGH CUT WAVENUMBER/2*PI CA CDB INTEGER I AMPLITUDE CUTOFF (DB) CA FCT REAL I FRACTION OF FX RANGE TO TAPER CA CA EJECT CAEND C*********************************************************************** C C LOCAL VARIABLES C C DPN - MINIMUM OF (P4-P3) AND (P2-P1) R*4 C F - ACTUAL FREQUENCY BEING CONSIDERED R*4 C FA,FB,FC,FD - USED IN DETERMINING FREQUENCY SEGMENTS R*4 C FXMAX- MAXIMUM FREQ*OFFSET COMPONENT R*4 C HMIN - THIS IS THE ACTUAL CUTOFF AMPLITUDE R*4 C IEF - ARRAY USED TO HOLD ENDING FREQ. INDEX FOR EACH C SEGMENT I*4 C IF - LOOP INDEX OVER THE FREQUENCIES I*4 C ISF - ARRAY USED TO HOLD STARTING FREQ. INDEX FOR EACH C SEGMENT I*4 C NSEG - NUMBER OF FREQUENCY SEGMMENTS TO USE (MAX OF 4) I*4 C PWR - LOG TO BASE 10 OF THE CUTOFF AMPLITUDE R*4 C P0 - AVERAGE RAY PARAMETER FOR TRAPEZOID R*4 C P1 - RAY PARAMETER FOR VELOCITY V1 R*4 C P2 - RAY PARAMETER FOR VELOCITY V2 R*4 C P3 - RAY PARAMETER FOR VELOCITY V3 R*4 C P4 - RAY PARAMETER FOR VELOCITY V4 R*4 C C*********************************************************************** C SUBROUTINE SAFF3DA( H,NXT,NF,DXT,XMX,SCF, * DF,F1,F2,F3,F4,V1,V2,V3,V4,SGN,ALB,CDB,FCT ) C IMPLICIT INTEGER (A-Z) C COMPLEX H(NXT,NF) REAL DXT(*),XMX(*),SCF(*) INTEGER ISF(4),IEF(4) C REAL ALA,ALB,AL0,AL1,AL2,AL3,AL4 REAL AN1,A1,A12,A2,A3,A34,A4 REAL B12,B34,C12,C34,DPN,DF,FCT REAL F,FA,FB,FC,FD,FE,FXMAX,F1,F2,F3,F4 REAL G12I,G12R,G34I,G34R,HI,HMIN,HR REAL PH3,PI,PWR,P0,P1,P2,P3,P4,SCT,SC1,S12,S34 REAL SAFF3DG,X,X1 C C --------------------------------------------------------------------- C PI = 4.0 * ATAN(1.0) C P1 = 1./V1 P2 = 1./V2 P3 = 1./V3 P4 = 1./V4 P0 = .5*(P1+P2-P3-P4) DPN = AMIN1( ABS(P4-P3),ABS(P2-P1) ) C C DETERMINE MAXIMUM FX VALUE C PWR = -CDB/20. HMIN = 10.**PWR PH3 = 16.*PI*ABS(P0)*DPN*DPN*HMIN PWR = -1./3. FXMAX = PH3**PWR C C --------------------------------------------------------------------- C C DETERMINE NUMBER OF SEGMENTS IN FREQUENCY AXIS C FA = F1 FB = V1*ALB ALA = .5*FB*(P1 + P2) ISF(1) = 1 C CINITIALIZE FOR POSSIBLE DIVIDE BY ZERO BY OPTIMIZER C FC = 1.0 FE = 1.0 FD = 0.0 CXX C C -------------------- C IF( FA.GT.FB ) THEN C IEF(1) = 0 NSEG = 1 ISF(2) = 1 C ELSE C IB = (FB-FA)/DF + 1 FB = FA + DF*(IB-1) C IF( FB.GE.F4 ) THEN C NSEG = 1 IEF(1) = NF GO TO 50 C ENDIF C IEF(1) = IB ISF(2) = IB + 1 C ENDIF C C -------------------- C ALA = .5*FB*(P1 + P2) FC = V2*ALA C IF( FA.GT.FC ) THEN C IEF(2) = 0 NSEG = 2 ISF(3) = 1 C ELSE C IC = (FC-FA)/DF + 1 FC = FA + DF*(IC-1) C IF( FC.GE.F4 ) THEN C NSEG = 2 IEF(2) = NF GO TO 50 C ENDIF C IEF(2) = IC ISF(3) = IC + 1 C ENDIF C C -------------------- C FD = V3*ALA ID = (FD-FA)/DF + 1 FD = FA + DF*(ID-1) C IF( FD.GE.F4 ) THEN C NSEG = 3 IEF(3) = NF GO TO 50 C ENDIF C C C -------------------- C IEF(3) = ID ISF(4) = ID + 1 FE = V4*ALB IE = (FE-FA)/DF + 1 FE = FA + DF*(IE-1) NSEG = 4 IEF(4) = NF C C -------------------------------------------------------------- C C FOR EACH SEGMENT COMPUTE THE TABLE OF FILTERS C 50 DO 1000 ISEG = 1,NSEG C DO 500 IF = ISF(ISEG),IEF(ISEG) F = F1 + DF*(IF-1) AL4 = F/V4 AL3 = F/V3 C IF( ISEG.EQ.1 ) THEN C AL1 = F/V1 AL2 = F/V2 XMX(IF) = FXMAX/F C ELSE IF( ISEG.EQ.2 ) THEN C AL1 = ALB AL2 = F/V2 XMX(IF) = FXMAX/F C ELSE IF( ISEG.EQ.3 ) THEN C AL1 = ALB AL2 = ALA XMX(IF) = FXMAX/FC C ELSE IF( ISEG.EQ.4 ) THEN C AL1 = ALB AL2 = ALA*(FE-F)/(FE-FD) + ALB*(F-FD)/(FE-FD) AL3 = AL2 XMX(IF) = FXMAX/FC C ENDIF C DXT(IF) = XMX(IF)/(NXT-1) C C OBTAIN THE TRAPEZOIDAL WEIGHTING FOR FREQUENCY F C SCF(IF) = SAFF3DG( F,F1,F2,F3,F4 ) C C NOW COMPUTE ENTRIES IN THE TABLE C AL0 = PI*( AL1 + AL2 - AL3 - AL4 ) H(1,IF) = CMPLX(AL0,0.) C DO 200 IX = 2,NXT X = (IX-1)*DXT(IF) A1 = 2.*AL1*X A2 = 2.*AL2*X A3 = 2.*AL3*X A4 = 2.*AL4*X A12 = A1 - A2 A34 = A3 - A4 C IF( A12.NE.1. ) THEN B12 = 1./(1.-A12*A12) C12 = COS(PI*A1) + COS(PI*A2) S12 = SIN(PI*A1) + SIN(PI*A2) ELSE B12 = .5 C12 = SIN(PI*A1) S12 = -COS(PI*A1) ENDIF C IF( A34.NE.1. ) THEN B34 = 1./(1.-A34*A34) C34 = COS(PI*A3) + COS(PI*A4) S34 = SIN(PI*A3) + SIN(PI*A4) ELSE B34 = .5 C34 = SIN(PI*A3) S34 = -COS(PI*A3) ENDIF C G12R = .5*B12*S12/X G12I = .5*B12*C12/X C G34R = .5*B34*S34/X G34I = .5*B34*C34/X C HR = (G12R - G34R) HI = (G12I - G34I) IF( SGN.EQ.-1 ) HI = -HI C 200 H(IX,IF) = CMPLX( HR,HI ) C C C APPLY TAPER C IF( FCT.GT.0. ) THEN C X1 = XMX(IF)*(1.-FCT)/DXT(IF) IX1 = X1 SC1 = PI/FLOAT(NXT-IX1) C DO 300 IX = IX1,NXT AN1 = FLOAT(IX-IX1) SCT = .5 + .5*COS(SC1*AN1) 300 H(IX,IF) = SCT*H(IX,IF) C ENDIF C 500 CONTINUE C 1000 CONTINUE C C --------------------------------------------------------------------- C RETURN END