CTITLESACNAY -- COHERENT NOISE ESTIMATION AND FILTERING 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR L. B. LIN 00020000 CA DESIGNER L. B. LIN 00030000 CA LANGUAGE FORTRAN 00040000 CA SYSTEM IBM AND CRAY 00050000 CA WRITTEN 07-20-85 00060000 C REVISED 06-05-86 LBL MODS FOR TWO MODE FILTERING 00070000 C REVISED 09-15-86 DJP ADDED THE BROADSIDE NOISE ATTENUATION 00080000 C OPTION 00090000 C REVISED 05-18-87 RDK IMPLEMENTED DJP'S CHANGES INTO 00100000 C PRODUCTION 00110000 C REVISED 03-08-89 LWC ADD SPATIAL WINDOW TAPER CAPABILITY 00120001 C REVISED 11-13-89 RDK FOR CRAY CFT77 COMPATIBILITY. 00121002 CA 00130000 CA 00140000 CA THIS ROUTINE HAS THREE ENTRY POINTS: 00150000 CA SACNAY FILTERS AND OUTPUTS A TRACE 00160000 CA SACNY2 STORES ATTRIBUTES OF AN INCOMING TRACE 00170000 CA SACNY3 CALCULATES COMPLEX NOISE FROM A PANEL OF TRACES 00180000 CA CALLING SEQUENCE IS DESCRIBED UNDER EACH ENTRY POINT. 00190000 CA 00200000 CA 00210000 C 00220000 C EJECT 00230000 C 00240000 SUBROUTINE SACNAY (RA, SA, OH, OTR, 00250000 * IPR, KPRTF, 00260000 * IXT, IXID, IXX, INOISE, IXK, IXST, IXET, 00270000 * LEN ,NF, IF1, M, N, NBL, LMZ, OSW, SR, 00280000 * TRC, CNTR, 00290000 * NCURVE, INOIS2, IXK2, BNAFLG) 00300000 C 00310000 C SACNAY FILTERS AND OUTPUTS A TRACE 00320000 C 00330000 C INPUT RA = RESERVED AREA OF BLANK COMMON FOR CNAY R4 00340000 C INPUT SA = SCRATCH AREA OF BLANK COMMON, WORKING SPACE R4 00350000 C OUTPUT OH = OUTPUT TRACE HEADER R4 00360000 C OUTPUT OTR = OUTPUT TRACE R4 00370000 C INPUT IPR = LOGICAL UNIT OF PRINTER I4 00380000 C INPUT KPRTF = KP AREA RETURN FLAG I4 00390000 C INPUT IXT = INDEX TO RA, POINTING TO TIME TRACE ARRAY I4 00400000 C INPUT IXID = INDEX TO RA, POINTING TO TRACE ID ARRAY I4 00410000 C INPUT IXX = INDEX TO RA, POINTING TO OFFSET ARRAY I4 00420000 C INPUT INOISE = INDEX TO RA, POINTING TO COMPLEX NOISE I4 00430000 C INPUT IXK = INDEX TO RA, POINTING TO WAVE NUMBER ARRAY I4 00440000 C INPUT IXST = INDEX TO RA, POINTING START TIME OF WINDOW I4 00450000 C INPUT IXET = INDEX TO RA, POINTING END TIME OF WINDOW I4 00460000 C INPUT LEN = LENGTH OF HEADER PLUS TIME TRACE I4 00470000 C INPUT NF = NUMBER OF FREQUENCIES I4 00480000 C INPUT IF1 = INDEX TO STARTING FREQUENCY I4 00490000 C INPUT M,N = FFT PARAMETERS: TRACE LENGTH N = 2**M I4 00500000 C INPUT NBL = NO. OF SAMPLES I4 00510000 C INPUT LMZ = LENGTH OF MERGING ZONE I4 00520000 C INPUT OSW = OUTPUT SWITCH R4 00530000 C INPUT SR = SAMPLING RATE R4 00540000 C INPUT TRC = INDEX TO THE TRACE TO BE FILTERED I4 00550000 C INPUT CNTR = INDEX TO THE CENTRAL TRACE I4 00560000 C INPUT BNAFLG = BROADSIDE NOISE ATTENUATION FLAG I4 00570000 C 00580000 C 00590000 DIMENSION RA(1) 00600000 DIMENSION SA(1) 00610000 DIMENSION OH(1) 00620000 DIMENSION OTR(1) 00630000 C 00640000 DIMENSION ISEG(1), IBEG(1), IEND(1) 00650000 INTEGER TICD 00690000 INTEGER TRC 00700000 INTEGER CNTR 00710000 INTEGER BNAFLG 00720000 C 00730000 REAL KX 00740000 REAL K1X 00750000 REAL K2X 00760000 REAL ATV 00770001 REAL SCF 00780000 REAL TAPWT 00790000 C 00800001 REAL TWAIT(55) 00810001 C 00820000 CALL ARMVE (RA(IXT+(TRC-1)*LEN), OH, LEN) 00830000 TRCID = RA(IXID+TRC-1) 00840000 C 00850000 IF (TRCID .EQ. 1.) THEN 00860000 C 00870000 X = RA(IXX+TRC-1) 00880000 XC = RA(IXX+CNTR-1) 00890000 DX = X - XC 00900000 CALL ARSET (SA(1), N+2, 0.) 00910000 IF = IF1 00920000 C 00930000 IF (DX .NE. 0.) GO TO 10 00940000 C 00950000 IF (NCURVE .LE. 1) THEN 00960000 CALL ARMVE (RA(INOISE), SA(IF1*2-1), 2*NF) 00970000 ELSE 00980000 CALL ARADF (RA(INOISE), RA(INOIS2), SA(IF1*2-1), 2*NF) 00990000 END IF 01000000 C 01010000 GO TO 30 01020000 C 01030000 CDIR$ IVDEP 01040000 10 DO 20 I = 1, NF 01050000 C 01060000 ARG = RA(IXK-1+I) * DX 01070000 COSARG = COS(ARG) 01080000 SINARG = SIN(ARG) 01090000 C 01100000 SA(IF*2-1) = RA(INOISE-1+I*2-1)*COSARG + 01110000 * RA(INOISE-1+I*2 )*SINARG 01120000 C 01130000 SA(IF*2 ) = RA(INOISE-1+I*2 )*COSARG - 01140000 * RA(INOISE-1+I*2-1)*SINARG 01150000 IF = IF + 1 01160000 C 01170000 20 CONTINUE 01180000 C 01190000 IF (NCURVE .LE. 1) GO TO 30 01200000 C 01210000 IF = IF1 01220000 CDIR$ IVDEP 01230000 DO 21 I = 1, NF 01240000 C 01250000 ARG = RA(IXK2-1+I) * DX 01260000 COSARG = COS(ARG) 01270000 SINARG = SIN(ARG) 01280000 C 01290000 SA(IF*2-1) = SA(IF*2-1) + 01300000 * RA(INOIS2-1+I*2-1)*COSARG + 01310000 * RA(INOIS2-1+I*2 )*SINARG 01320000 C 01330000 SA(IF*2 ) = SA(IF*2 ) + 01340000 * RA(INOIS2-1+I*2 )*COSARG - 01350000 * RA(INOIS2-1+I*2-1)*SINARG 01360000 IF = IF + 1 01370000 C 01380000 21 CONTINUE 01390000 C 01400000 30 CALL S2DFI2 (M, SA(1), * 40 )01410000 GO TO 50 01420000 C 01430000 40 KPRTF = -1 01440000 WRITE (IPR, 9010) 01450000 RETURN 01460000 C 01470000 C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 01480000 C 01490000 50 ISB = INT(RA(IXST+TRC-1) / SR + 1.001) 01500002 LWB = INT((RA(IXET+TRC-1) - RA(IXST+TRC-1)) / SR + 0.001) 01510002 SC = 0. 01520000 DS = 1./LMZ 01530000 J = ISB 01540000 C 01550000 DO 60 I = 1, LMZ 01560000 SA(J) = SC*SA(J) 01570000 SC = SC + DS 01580000 J = J + 1 01590000 C 01600000 60 CONTINUE 01610000 C 01620000 J = ISB + LWB - LMZ 01630000 SC = 1. 01640000 C 01650000 DO 70 I = 1, LMZ 01660000 SA(J) = SC*SA(J) 01670000 SC = SC - DS 01680000 J = J + 1 01690000 C 01700000 70 CONTINUE 01710000 C 01720000 IF (OSW .NE. 1.) THEN 01730000 CALL ARSBF (OTR(ISB), SA(ISB), OTR(ISB), LWB) 01740000 C 01750000 ELSE 01760000 CALL ARSET (OTR, NBL, 0.) 01770000 CALL ARMVE (SA(ISB), OTR(ISB), LWB) 01780000 END IF 01790000 C 01800000 C REVERSE THE SCALING DONE IN THE BNA OPTION 01810000 C 01820000 IF (BNAFLG .NE. 99) THEN 01830000 SCR = 1.0 / SQRT(RA(IXX+TRC-1)) 01840000 C 01850000 DO 80 I = 1, NBL 01860000 OTR(I) = OTR(I) * SCR 01870000 C 01880000 80 CONTINUE 01890000 C 01900000 END IF 01910000 C 01920000 C ********************************************************* 01930000 C 01940000 ELSE 01950000 C 01960000 CALL ARSET (OTR, NBL, 0.0) 01970000 C 01980000 ENDIF 01990000 C 02000000 RETURN 02010000 C 02020000 C 02030000 C 02040000 ENTRY SACNY2 (RA, SA, OH, OTR, 02050000 * IPR, KPRTF, 02060000 * IXOFST, IXSTM, IXETM, 02070000 * IXX, IXID, IXST, IXET, IXD, IXT, 02080000 * XDST, TICD, NXSE, SR, M, N, IF1, NF, LEN, 02090000 * KTRC, BNAFLG, BNXC, BNYC) 02100000 C 02110000 C SACNY2 STORES ATTRIBUTES OF AN INCOMING TRACE 02120000 C 02130000 C INPUT RA = BLANK COMMON RESERVED FOR CNAY R4 02140000 C INPUT SA = SCRATCH AREA OF BLANK COMMON, WORKING SPACE R4 02150000 C INPUT OH = INCOMING TRACE HEADER R4 02160000 C INPUT OTR = INCOMING TRACE R4 02170000 C INPUT IPR = LOGICAL UNIT OF PRINTER I4 02180000 C INPUT KPRTF = KP AREA RETURN FLAG I4 02190000 C INPUT IXOFST = INDEX TO RA, POINTING TO X-ST-ET TRIPLET I4 02200000 C INPUT IXSTM = INDEX TO RA, POINTING TO X-ST-ET TRIPLET I4 02210000 C INPUT IXETM = INDEX TO RA, POINTING TO X-ST-ET TRIPLET I4 02220000 C INPUT IXX = INDEX TO RA, POINTING TO OFFSET ARRAY I4 02230000 C INPUT IXID = INDEX TO RA, POINTING TO TRACE ID ARRAY I4 02240000 C INPUT IXST = INDEX TO RA, POINTING START TIME OF WINDOW I4 02250000 C INPUT IXET = INDEX TO RA, POINTING END TIME OF WINDOW I4 02260000 C INPUT IXD = INDEX TO RA, POINTING TO FREQ TRACE ARRAY I4 02270000 C INPUT IXT = INDEX TO RA, POINTING TO TIME TRACE ARRAY I4 02280000 C INPUT XDST = OFFSET DISTANCE R4 02290000 C INPUT TICD = TRACE ID I4 02300000 C INPUT NXSE = NUMBER OF X-ST-ET TRIPLETS I4 02310000 C INPUT SR = SAMPLING RATE R4 02320000 C INPUT M,N = FFT PARAMETERS: TRACE LENGTH N = 2**M I4 02330000 C INPUT IF1 = INDEX TO STARTING FREQUENCY I4 02340000 C INPUT NF = NUMBER OF FREQUENCIES I4 02350000 C INPUT LEN = LENGTH OF HEADER PLUS TIME TRACE I4 02360000 C INPUT KTRC = INDEX OF TRACE ACCUMULATION I4 02370000 C INPUT BNAFLG = BROADSIDE NOISE ATTENUATION FLAG I4 02380000 C INPUT BNXC = X-COORDINATE OF THE SECONDARY NOISE SOURCE R4 02390000 C INPUT BNYC = Y-COORDINATE OF THE SECONDARY NOISE SOURCE R4 02400000 C OUTPUT NUMEROUS ATTRIBUTES OF TRACE STORED IN RA 02410000 C 02420000 C 02430000 RA(IXX+KTRC-1) = XDST 02440000 RA(IXID+KTRC-1) = TICD 02450000 C 02460000 C COMPUTE THE EFFECTIVE OFFSET FOR THE BNA OPTION 02470000 C 02480000 IF (BNAFLG .NE. 99) THEN 02490000 CALL USRTHV (OH, 'THRXC ', IRXC) 02500000 CALL USRTHV (OH, 'THRYC ', IRYC) 02510000 RXC = IRXC 02520000 RYC = IRYC 02530000 DRX = RXC - BNXC 02540000 DRY = RYC - BNYC 02550000 C 02560000 RA(IXX+KTRC-1) = SQRT (DRX * DRX + DRY * DRY) 02570000 SCR = SQRT (RA(IXX+KTRC-1)) 02580000 C 02590000 C SCALE TRACE BY SQRT OF THE NEW OFFSET 02600000 C 02610000 DO 90 I = 1, NBL 02620000 OTR(I) = OTR(I) * SCR 02630000 C 02640000 90 CONTINUE 02650000 C 02660000 END IF 02670000 C 02680000 XA = RA(IXX+KTRC-1) 02690000 CALL MTRPLT (NXSE, RA(IXOFST), RA(IXSTM), RA(IXETM), 02700000 * XA, STMA, ETMA) 02710000 RA(IXST+KTRC-1) = STMA 02720000 RA(IXET+KTRC-1) = ETMA 02730000 C 02740000 ISB = INT(STMA / SR + 1.001) 02750002 LWB = INT((ETMA - STMA) / SR + 0.001) 02760002 C 02770000 IF (TICD .EQ. 1) THEN 02780000 CALL ARSET (SA(1), N, 0.) 02790000 CALL ARMVE (OTR(ISB), SA(ISB), LWB) 02800000 CALL S2DFT2 (M, SA(1), * 100 )02810000 GO TO 110 02820000 C 02830000 100 KPRTF = -1 02840000 WRITE (IPR, 9000) 02850000 RETURN 02860000 C 02870000 110 IS = 2*IF1 - 1 02880000 CALL ARMVE (SA(IS), RA(IXD+(KTRC-1)*2*NF), 2*NF) 02890000 C 02900000 C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 02910000 ELSE 02920000 CALL ARSET (RA(IXD+(KTRC-1)*2*NF), 2*NF, 0.0) 02930000 ENDIF 02940000 C 02950000 CALL ARMVE (OH, RA(IXT+(KTRC-1)*LEN), LEN) 02960000 C 02970000 RETURN 02980000 C 02990000 C 03000000 C 03010000 C 03020000 ENTRY SACNY3 (RA, SA, 03030000 * IXX, IXID, IXD, IXK, INOISE, 03040000 * N2U, NF, MTPER, 03050000 * CNTR, NCURVE, 03060000 * MSG, ISEG, IBEG, IEND, INOIS2, IXK2, ATV) 03070001 C 03080000 C SACNY3 CALCULATES COMPLEX NOISE FROM A PANEL OF TRACES 03090000 C 03100000 C INPUT RA = BLANK COMMON RESERVED FOR CNAY R4 03110000 C INPUT SA = SCRATCH AREA R4 03120000 C INPUT IXX = INDEX TO RA, POINTING TO OFFSET ARRAY I4 03130000 C INPUT IXID = INDEX TO RA, POINTING TO TRACE ID ARRAY I4 03140000 C INPUT IXD = INDEX TO RA, POINTING TO FREQ TRACE ARRAY I4 03150000 C INPUT IXK = INDEX TO RA, POINTING TO WAVE NUMBER ARRAY I4 03160000 C INPUT INOISE = INDEX TO RA, POINTING TO COMPLEX NOISE I4 03170000 C INPUT N2U = NUMBER OF TRACES IN THE COMPUTATION PANEL I4 03180000 C INPUT NF = NUMBER OF FREQUENCIES I4 03190000 C INPUT MTPER = NUMBER OF FREQUENCIES TO APPLY TAPER I4 03200000 C INPUT CNTR = POINTER TO THE CENTRAL TRACE I4 03210000 C INPUT ATV = PARAMETER USED IN EQUATION TO FIND WEIGHT R4 03220001 C OUTPUT RA(INOISE) CONTAINS COMPLEX NOISE 03230000 C 03240000 C 03250000 XC = RA(IXX+CNTR-1) 03260000 M2U = 0 03270000 C 03280001 IF (NCURVE .GE. 2) THEN 03290001 DO 115 J = 1,N2U 03300001 IF (RA(IXID-1+J) .EQ. 1.) M2U = M2U + 1 03310001 115 CONTINUE 03320001 IF (M2U .EQ. 0) RETURN 03330001 C 03340001 SCF = 1.0 / M2U 03350001 C 03360001 GO TO 200 03370001 ENDIF 03380001 C 03390001 STAPWT = 0.0 03400001 N2UH = (N2U + 1) / 2 03410001 NXDIF = CNTR - N2UH 03420001 C 03430000 DO 120 J = 1,N2U 03440001 IF (RA(IXID-1+J) .EQ. 1.) THEN 03450001 M2U = M2U + 1 03460001 IF (ATV .EQ. 0.0) THEN 03470001 TWAIT(J) = 1.0 03480001 STAPWT = STAPWT + TWAIT(J) 03490001 GO TO 120 03500001 ENDIF 03510001 JJ = J 03520001 IF (J .GT. N2UH) THEN 03530001 JJ = N2U - J + 1 03540001 ENDIF 03550001 TWAIT(J) = EXP(-ATV * (JJ-N2UH)**2) 03560001 STAPWT = STAPWT + TWAIT(J) 03570001 C 03580000 ELSE 03590001 TWAIT(J) = 0.0 03600001 ENDIF 03610001 120 CONTINUE 03620001 C 03630001 C 03800000 IF (M2U .EQ. 0) RETURN 03810001 C 03820000 C 03830000 C CALCULATION OF NOISE FOR ONE MODE CASE 03840000 C 03850000 CDIR$ IVDEP 03860000 DO 130 I = 1, NF 03870000 SA( I) = 0.0 03880000 SA(NF+I) = 0.0 03890000 C 03900000 130 CONTINUE 03910000 C 03920000 C 03930000 DO 150 J = 1, N2U 03940000 C 03950000 JX = J - NXDIF 03960001 IF (JX .LE. 0) JX = JX + N2U 03970001 IF (JX .GT. N2U) JX = JX - N2U 03980001 XMXC = RA(IXX-1+J) - XC 04010000 JM1NF = (J-1) * NF 04020000 C 04030000 CDIR$ IVDEP 04040000 DO 140 I = 1, NF 04050000 C 04060000 ARG = RA(IXK-1+I) * XMXC 04070000 COSARG = COS(ARG) 04080000 SINARG = SIN(ARG) 04090000 C 04100000 IA = I + JM1NF 04110000 IA2 = IA * 2 04120000 IA2M1 = IA2 - 1 04130000 AAR = RA(IXD-1+IA2M1) 04140000 AAI = RA(IXD-1+ IA2) 04150000 C 04160000 DUMR = AAR * COSARG - AAI * SINARG 04170000 DUMI = AAR * SINARG + AAI * COSARG 04180000 C 04190000 SA( I) = SA( I) + DUMR * TWAIT(JX) 04200001 SA(NF+I) = SA(NF+I) + DUMI * TWAIT(JX) 04210001 C 04220000 140 CONTINUE 04230000 C 04240000 150 CONTINUE 04250000 SCF = 1.0 / STAPWT 04260000 C 04270001 DO 160 I = 1, NF 04330000 I2M1 = I * 2 - 1 04340000 RA(INOISE-1+I2M1) = SA( I) * SCF 04350000 C 04360000 160 CONTINUE 04370000 C 04380000 DO 170 I = 1, NF 04390000 I2 = I * 2 04400000 RA(INOISE-1+I2 ) = SA(NF+I) * SCF 04410000 C 04420000 170 CONTINUE 04430000 C 04440001 C 04520000 C***************************************************************** 04530000 C 04540000 C APPLY FREQUENCY TAPER 04550000 C 04560000 IF (MTPER .GT. 0) THEN 04570000 SC = 0. 04580000 DS = 1. / MTPER 04590000 CDIR$ IVDEP 04600000 DO 180 I = 1, MTPER 04610000 I2 = I * 2 04620000 RA(INOISE-1+I2-1) = SC * RA(INOISE-1+I2-1) 04630000 RA(INOISE-1+I2 ) = SC * RA(INOISE-1+I2 ) 04640000 SC = SC + DS 04650000 C 04660000 180 CONTINUE 04670000 C 04680000 SC = 1. 04690000 CDIR$ IVDEP 04700000 DO 190 I = 1, MTPER 04710000 I2 = I * 2 04720000 SC = SC - DS 04730000 RA(INOISE+2*NF-2*MTPER-1+I2-1) = 04740000 * SC*RA(INOISE+2*NF-2*MTPER-1+I2-1) 04750000 C 04760000 RA(INOISE+2*NF-2*MTPER-1+I2 ) = 04770000 * SC*RA(INOISE+2*NF-2*MTPER-1+I2 ) 04780000 C 04790000 190 CONTINUE 04800000 C 04810000 ENDIF 04820000 C 04830000 RETURN 04840000 C 04850000 C CALCULATION OF NOISE FOR TWO MODE CASE 04860000 C 04870000 200 CALL ARSET (RA(INOISE), 2*NF, 0.) 04880000 CALL ARSET (RA(INOIS2), 2*NF, 0.) 04890000 C 04900000 DO 450 I = 1, MSG 04910000 C 04920000 IS = IBEG(I) 04930000 IE = IEND(I) 04940000 C 04950000 IF (ISEG(I).EQ.2 .OR. ISEG(I).EQ.3) THEN 04960000 C 04970000 C CALL SACNA1 (RA(IXD),RA(IXX),N2U,NF,SCF,XC,IBEG(I),IEND(I), 04980000 C * RA(IXK), RA(INOISE), SA(1),SA(1+NF) ) 04990000 C ELSE IF (ISEG(I) .EQ. 3) THEN 05000000 C CALL SACNA1 (RA(IXD),RA(IXX),N2U,NF,SCF,XC,IBEG(I),IEND(I), 05010000 C * RA(IXK2), RA(INOIS2), SA(1),SA(1+NF)) 05020000 C 05030000 IF (ISEG(I).EQ.2) IXKC = IXK 05040000 IF (ISEG(I).EQ.2) INOISC = INOISE 05050000 IF (ISEG(I).EQ.3) IXKC = IXK2 05060000 IF (ISEG(I).EQ.3) INOISC = INOIS2 05070000 C 05080000 CDIR$ IVDEP 05090000 DO 220 J = IS, IE 05100000 SA(J) = 0. 05110000 220 SA(NF+J) = 0. 05120000 C 05130000 DO 240 J = 1, N2U 05140000 X = RA(IXX+J-1) - XC 05150000 C 05160000 CDIR$ IVDEP 05170000 DO 230 K = IS, IE 05180000 KX = RA(IXKC+K-1)*X 05190000 C 05200000 IDXI = ((J-1)*NF + K)*2 05210000 IDXR = IDXI - 1 05220000 DR = RA(IXD+IDXR-1) 05230000 DI = RA(IXD+IDXI-1) 05240000 C 05250000 COSKX = COS(KX) 05260000 SINKX = SIN(KX) 05270000 C 05280000 SA(K) = SA(K) + DR*COSKX - DI*SINKX 05290000 SA(NF+K) = SA(NF+K) + DR*SINKX + DI*COSKX 05300000 230 CONTINUE 05310000 C 05320000 240 CONTINUE 05330000 C 05340000 CDIR$ IVDEP 05350000 DO 250 J = IS, IE 05360000 RA(INOISC+J*2-2) = SA( J)*SCF 05370000 RA(INOISC+J*2-1) = SA(NF+J)*SCF 05380000 250 CONTINUE 05390000 C 05400000 ELSE IF (ISEG(I) .EQ. 4) THEN 05410000 C 05420000 C CALL SACNA2 (RA(IXD),RA(IXX),N2U,NF,SCF,XC,IBEG(I),IEND(I), 05430000 C * RA(IXK), RA(IXK2),RA(INOISE),RA(INOIS2), 05440000 C * SA(1+0*NF),SA(1+1*NF),SA(1+2*NF), 05450000 C * SA(1+3*NF),SA(1+4*NF),SA(1+5*NF)) 05460000 C 05470000 C 05480000 CDIR$ IVDEP 05490000 DO 260 J = IS, IE 05500000 SA(0*NF+J) = 0. 05510000 SA(1*NF+J) = 0. 05520000 SA(2*NF+J) = 0. 05530000 SA(3*NF+J) = 0. 05540000 SA(4*NF+J) = 0. 05550000 SA(5*NF+J) = 0. 05560000 260 CONTINUE 05570000 C 05580000 DO 300 J = 1, N2U 05590000 X = RA(IXX+J-1) - XC 05600000 C 05610000 CDIR$ IVDEP 05620000 DO 280 K = IS, IE 05630000 K1X = RA(IXK +K-1)*X 05640000 K2X = RA(IXK2+K-1)*X 05650000 IDXI= ((J-1)*NF + K)*2 05660000 IDXR= IDXI - 1 05670000 DR = RA(IXD+IDXR-1) 05680000 DI = RA(IXD+IDXI-1) 05690000 C 05700000 COSK1X = COS(K1X) 05710000 COSK2X = COS(K2X) 05720000 SINK1X = SIN(K1X) 05730000 SINK2X = SIN(K2X) 05740000 C 05750000 SA(0*NF+K) = SA(0*NF+K) + DR*COSK1X - DI*SINK1X 05760000 SA(1*NF+K) = SA(1*NF+K) + DR*SINK1X + DI*COSK1X 05770000 SA(2*NF+K) = SA(2*NF+K) + DR*COSK2X - DI*SINK2X 05780000 SA(3*NF+K) = SA(3*NF+K) + DR*SINK2X + DI*COSK2X 05790000 C 05800000 SA(4*NF+K) = SA(4*NF+K) + COSK1X*COSK2X + SINK1X*SINK2X 05810000 SA(5*NF+K) = SA(5*NF+K) + SINK1X*COSK2X - COSK1X*SINK2X 05820000 C 05830000 280 CONTINUE 05840000 C 05850000 300 CONTINUE 05860000 C 05870000 CDIR$ IVDEP 05880000 DO 400 J = IS, IE 05890000 SA(0*NF+J) = SA(0*NF+J)*SCF 05900000 SA(1*NF+J) = SA(1*NF+J)*SCF 05910000 SA(2*NF+J) = SA(2*NF+J)*SCF 05920000 SA(3*NF+J) = SA(3*NF+J)*SCF 05930000 SA(4*NF+J) = SA(4*NF+J)*SCF 05940000 SA(5*NF+J) = SA(5*NF+J)*SCF 05950000 C 05960000 FAC = 1./(1. - (SA(4*NF+J)**2 + SA(5*NF+J)**2)) 05970000 C 05980000 RA(INOISE+J*2-2) = (SA(0*NF+J)-SA(4*NF+J)*SA(2*NF+J)+ 05990000 * SA(5*NF+J)*SA(3*NF+J))*FAC 06000000 RA(INOISE+J*2-1) = (SA(1*NF+J)-SA(4*NF+J)*SA(3*NF+J)- 06010000 * SA(5*NF+J)*SA(2*NF+J))*FAC 06020000 RA(INOIS2+J*2-2) = (SA(2*NF+J)-SA(4*NF+J)*SA(0*NF+J)- 06030000 * SA(5*NF+J)*SA(1*NF+J))*FAC 06040000 RA(INOIS2+J*2-1) = (SA(3*NF+J)-SA(4*NF+J)*SA(1*NF+J)+ 06050000 * SA(5*NF+J)*SA(0*NF+J))*FAC 06060000 400 CONTINUE 06070000 C 06080000 END IF 06090000 C 06100000 450 CONTINUE 06110000 C 06120000 C***************************************************************** 06130000 C 06140000 C APPLY FREQUENCY TAPER TO BOTH NOISE MODES 06150000 C 06160000 IF (MTPER .GT. 0) THEN 06170000 SC = 0. 06180000 DS = 1. / MTPER 06190000 CDIR$ IVDEP 06200000 DO 500 I = 1, MTPER 06210000 I2 = I * 2 06220000 RA(INOISE-1+I2-1) = SC * RA(INOISE-1+I2-1) 06230000 RA(INOISE-1+I2 ) = SC * RA(INOISE-1+I2 ) 06240000 RA(INOIS2-1+I2-1) = SC * RA(INOIS2-1+I2-1) 06250000 RA(INOIS2-1+I2 ) = SC * RA(INOIS2-1+I2 ) 06260000 SC = SC + DS 06270000 C 06280000 500 CONTINUE 06290000 C 06300000 SC = 1. 06310000 CDIR$ IVDEP 06320000 DO 600 I = 1, MTPER 06330000 I2 = I * 2 06340000 SC = SC - DS 06350000 RA(INOISE+2*NF-2*MTPER-1+I2-1) = 06360000 * SC*RA(INOISE+2*NF-2*MTPER-1+I2-1) 06370000 C 06380000 RA(INOISE+2*NF-2*MTPER-1+I2 ) = 06390000 * SC*RA(INOISE+2*NF-2*MTPER-1+I2 ) 06400000 C 06410000 RA(INOIS2+2*NF-2*MTPER-1+I2-1) = 06420000 * SC*RA(INOIS2+2*NF-2*MTPER-1+I2-1) 06430000 C 06440000 RA(INOIS2+2*NF-2*MTPER-1+I2 ) = 06450000 * SC*RA(INOIS2+2*NF-2*MTPER-1+I2 ) 06460000 C 06470000 600 CONTINUE 06480000 C 06490000 ENDIF 06500000 C 06510000 RETURN 06520000 C 06530000 C FORMAT STATEMENTS 06540000 C 06550000 9000 FORMAT(1X,'ERROR CALLING "S2DFT2" IN SACNAY') 06560000 C 06570000 9010 FORMAT(1X,'ERROR CALLING "S2DFI2" IN SACNY2') 06580000 C 06590000 END 06600000