CTITLEARSQM -- SLIDING SUMS OF SQUARES 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR STU NELAN 00020000 CA DESIGNER BILL CACKLER (FROM SPARC ON SEL) 00030000 CA LANGUAGE FORTRAN 00040000 CA SYSTEM IBM AND CRAY 00050001 CA WRITTEN 01-15-86 00060000 C REVISED 07-11-86 ESN. FOR USE ON THE IBM. 00070001 CA 00080000 CA 00090000 CA CALL ARSQM (A, B, N, WL) 00100000 CA INPUT A = REAL ARRAY R4 00110000 CA OUTPUT B = REAL ARRAY R4 00120000 CA INPUT N = NUMBER ELEMENTS IN A I4 00130000 CA INPUT WL = WINDOW LENGTH IN POINTS I4 00140000 CA 00150000 CA 00160000 CA ARSQM SLIDES A WINDOW OF 'WL' POINTS ALONG THE INPUT ARRAY A 00170000 CA AND WRITES IN ARRAY B THE SUM OF SQUARED VALUES DIVIDED BY 00180000 CA THE NUMBER OF VALUES IN THE WINDOW, NOT COUNTING TWO OR MORE 00190000 CA ZEROS. THE OUTPUT LENGTH OF B WILL BE EQUAL TO N - WL + 1. 00200000 CAEND 00210000 SUBROUTINE ARSQM (A, B, N, WL) 00220000 C 00230000 REAL A (1) 00240000 REAL B (1) 00250000 C 00260000 INTEGER WINCNT 00270000 INTEGER WL 00280000 INTEGER ZEROCT 00290000 C 00300000 C INITIALIZATION 00310000 C 00320000 IF (WL .LT. 3 .OR. WL. GT. N) GO TO 300 00330000 WINCNT = WL 00340000 ZEROCT = 0 00350000 ASUM = 0. 00360000 SUM = 0. 00370000 C 00380000 C COMPUTE THE SUM OF THE SQUARES AND THE 00390000 C NUMBER OF POINTS TO USE IN FINDING THE 00400000 C AVERAGE OVER THE FIRST WINDOW. 00410000 C 00420000 DO 100 I = 1, WL 00430000 SUM = SUM + A(I) * A(I) 00440000 IF (A(I) .EQ. 0.0) ZEROCT = ZEROCT + 1 00450000 IF (A(I) .NE. 0.0 .AND. ZEROCT .GT. 1) WINCNT = WINCNT - ZEROCT 00460000 IF (A(I) .NE. 0.0) ZEROCT = 0 00470000 100 CONTINUE 00480000 C 00490000 C IF WINDOW ENDED WITH 2 OR MORE ZEROS, 00500000 C ADJUST FOR THEM AND THEN COMPUTE THE AVERAGE 00510000 C SQUARE FOR THE WINDOW AND PUT IN PROPER 00520000 C PLACE IN THE OUTPUT ARRAY. 00530000 C 00540000 IF (ZEROCT .GT. 1) WINCNT = WINCNT - ZEROCT 00550000 IF (SUM .NE. 0.0) ASUM = SUM / WINCNT 00560000 B(1) = ASUM 00570000 C 00580000 C IF ONLY 1 WINDOW TO DO, RETURN 00590000 C 00600000 IF (N .EQ. WL) GO TO 300 00610000 ISTART = WL + 1 00620000 C 00630000 C COMPUTE AVERAGE SQUARE OVER EACH REMAINING 00640000 C WINDOW BY SIMPLY ADDING THE NEXT SQUARE AND 00650000 C SUBTRACTING THE FIRST SQUARE OF THE PREVIOUS 00660000 C WINDOW. THE NUMBER OF POINTS TO DIVIDE BY 00670000 C MUST BE ADJUSTED BY LOOKING AT THE LAST THREE 00680000 C POINTS OF THE NEW WINDOW AND THE FIRST THREE 00690000 C POINTS OF THE PREVIOUS WINDOW. 00700000 C 00710000 DO 200 I = ISTART, N 00720000 IADD = 0 00730000 ISUB = 0 00740000 ASUM = 0. 00750000 SUM = SUM + A(I) * A(I) - A(I-WL) * A(I-WL) 00760000 IF (A(I-WL) .EQ. 0.0 .AND. A(I-WL+1) .EQ. 0.0) IADD = 1 00770000 IF (A(I) .EQ. 0.0 .AND. A(I-1) .EQ. 0.0) ISUB = 1 00780000 IF (IADD .EQ. 1 .AND. A(I-WL+2) .NE. 0.0) IADD = 2 00790000 IF (ISUB .EQ. 1 .AND. A(I-2) .NE. 0.0) ISUB = 2 00800000 WINCNT = WINCNT + IADD - ISUB 00810000 IF (SUM .NE. 0.0) ASUM = SUM / WINCNT 00820000 B(I-WL+1) = ASUM 00830000 200 CONTINUE 00840000 C 00850000 300 CONTINUE 00860000 C 00870000 RETURN 00880000 C 00890000 END 00900000