*PROCESS; 00010000 LSPYSMO: 00020000 PROC ( F, DX, NUMBER_OF_SAMPLES, NPTS_PER_FIT, NORDER,NPTS_PER_JUMP);00030000 % SKIP(001); 00040000 /********************************************************************00050000 * *00060000 * THIS ROUTINE SMOOTHES THE INPUT FUNCTION F BY FITTING A LEAST *00070000 * SQUARES POLYNOMIAL THROUGH SUCCESSIVE SETS OF POINTS OF F AND *00080000 * RESAMPLES F AT THESE POINTS USING THE POLYNOMIALS FOUND. *00090000 * *00100000 * INPUT PARAMETERS: *00110001 * *00120001 * F(*) THE FLOATING POINT ARRAY OF VALUES TO BE *00130001 * SMOOTHED. *00140001 * DX THE SAMPLE INTERVAL BETWEEN THE SAMPLES *00150001 * IN F. *00160001 * NUMBER_OF_SAMPLES THE NUMBER OF SAMPLES IN THE ARRAY F. *00170001 * NPTS_PER_FIT THE NUMBER OF SAMPLE POINTS TO BE USED *00180001 * AT ONE TIME TO FIT F. *00190001 * NORDER THE ORDER OF THE POLYNOMIAL TO FIT F. *00200001 * NPTS_PER_JUMP THE NUMER OF SAMPLE POINTS TO MOVE FORWARD*00210001 * FOR THE NEXT POLYNOMIAL FIT OF F. *00220001 * *00230001 * OUTPUT: *00240001 * *00250001 * F(*) THE SMOOTHED VALUES OF F ARE RETURNED *00260001 * IN F. *00270001 * BY: *00280001 * L. DAVID JONES *00290001 * SEPT. 1977 *00300001 ********************************************************************/00310001 % SKIP(001); 00320000 DCL F(*) FLOAT BIN; 00330001 DCL (A(*,*), B(*,*)) CONTROLLED FLOAT BIN; 00340000 DCL X FLOAT BIN; 00350000 DCL ( NUMBER_OF_SAMPLES, NPTS_PER_FIT, NORDER, NPTS_PER_JUMP, 00360001 IPTS_PER_FIT, NTERMS, NUM1, NTERMSM2, IS 00370000 ) FIXED BIN(31); 00380000 DECLARE INDEX FIXED BIN; 00390000 % SKIP(001); 00400000 /* FIT THE FIRST SERIES OF POINTS OF THE VECTOR 00410000 ARRAY F TO AN NTH ORDER CURVE (TWO EXCESS POINTS 00420000 ARE INCLUDED IN THIS FIT). 00430000 */00440000 % SKIP(001); 00450000 ALLOCATE A(NPTS_PER_FIT+2,NORDER+1), B(NPTS_PER_FIT+2,1); 00460000 NUM1 = 1; 00470000 IPTS_PER_FIT = NPTS_PER_FIT + 2; 00480000 NTERMS = NORDER + 1; 00490000 NTERMSM2 = NTERMS - 2; 00500000 INDEX = 1; 00510000 CALL POINT_LSTSQRS_FITTER; 00520000 DO I=1 TO 4; 00530000 X = (I-1)*DX; 00540000 DUM = B(1,1); 00550000 DO J=2 TO NTERMS; 00560000 DUM = DUM + B(J,1)*X**(J-1); 00570000 END; 00580000 F(I) = DUM; 00590000 END; 00600000 F2 = F(4); 00610000 % SKIP(001); 00620000 /* GO THROUGH REST OF F FITTING NPTS_PER_FIT AT A TIME 00630000 TO AN NORDER(TH) ORDER CURVE UP TO JUST BEFORE LAST 00640000 SET OF POINTS. THIS LAST SET OF POINTS IS HANDLED 00650000 AS A SPECIAL CASE SO THAT ONE DOESN'T RUN OFF THE 00660000 END OF THE VECTOR ARRAY F. 00670000 */00680000 % SKIP(001); 00690000 IPTS_PER_FIT = NPTS_PER_FIT; 00700000 DO INDEX=3 TO NUMBER_OF_SAMPLES-NPTS_PER_FIT+1 BY NPTS_PER_JUMP; 00710000 CALL SPLINE_LSTSQRS_FITTER; 00720000 DO I=INDEX+2 TO INDEX + NPTS_PER_JUMP +1; 00730000 X = (I-INDEX )*DX; 00740000 DUM = B(1,1); 00750000 DO J=2 TO NTERMS; 00760000 DUM = DUM + B(J,1)*X**(J-1); 00770000 END; 00780000 F(I) = DUM; 00790000 END; 00800000 F2 = DUM; 00810000 END; 00820000 I = IS; % SKIP(001); 00840000 /* LAST SET OF POINTS TREATED AS A SPECIAL CASE 00850000 */00860000 % SKIP(001); 00870000 IF INDEX - NPTS_PER_JUMP ^= NUMBER_OF_SAMPLES - NPTS_PER_FIT + 1 THEN DO INDEX = NUMBER_OF_SAMPLES - NPTS_PER_FIT + 1; CALL SPLINE_LSTSQRS_FITTER; DO I = IS TO NUMBER_OF_SAMPLES; X = (I-IS)*DX; DUM = B(1,1); DO J = 2 TO NTERMS; DUM = DUM + B(J,1)*X**(J-1); END; F(I) = DUM; END; END; ELSE DO I=IS TO NUMBER_OF_SAMPLES; 00880000 X = (I-IS) * DX; DUM = B(1,1); 00900000 DO J=2 TO NTERMS; 00910000 DUM = DUM + B(J,1)*X**(J-1); 00920000 END; 00930000 F(I) = DUM; 00940000 END; 00950000 % PAGE; 00960000 POINT_LSTSQRS_FITTER: PROC; 00970000 % SKIP(001); 00980000 /****************************************************************00990000 * *01000000 * THIS ROUTINE FITS A LEAST SQUARES POLYNOMIAL OF ORDER *01010000 * NORDER, THROUGH 'IPTS_PER_FIT' NUMBER OF POINTS. IT THEN *01020000 * SAVES F, X, AND F' AT THE LAST TWO OF THOSE POINTS FOR *01030000 * SUBSEQUENT USE IN FITS USING THE PROC SPLINE_LSTSQRS_FITTER *01040000 * *01050000 ****************************************************************/01060000 % SKIP(001); 01070000 DO I=1 TO IPTS_PER_FIT; 01080000 X = (I-1)*DX; 01090000 B(I,1) = F(I); 01100000 A(I,1) = 1.E0; 01110000 DO J=2 TO NTERMS; 01120000 A(I,J) = X**(J-1); 01130000 END; 01140000 END; 01150000 CALL MLSQ ( A, B, IPTS_PER_FIT, NTERMS, NUM1 ); 01160000 XD1 = 2*DX; 01170000 XD2 = 3*DX; 01180000 F1 = F(3); 01190000 F2S = F(4); 01200000 FP1 = 0.0; 01210000 FP2 = 0.0; 01220000 DO J=2 TO NTERMS; 01230000 FP1 = FP1 + (J-1)*B(J,1)*XD1**(J-2); 01240000 FP2 = FP2 + (J-1)*B(J,1)*XD2**(J-2); 01250000 END; 01260000 END POINT_LSTSQRS_FITTER; 01270000 % PAGE; 01280000 SPLINE_LSTSQRS_FITTER: PROC; 01290000 % SKIP(001); 01300000 /****************************************************************01310000 * *01320000 * THIS ROUTINE FITS A LEAST SQUARES POLYNOMIAL THROUGH *01330000 * 'IPTS_PER_FIT' NUMBER OF POINTS PLUS A SLOPE CONDITION *01340000 * FROM THE LAST FIT AT THE FIRST OF THE CURRENT GROUP OF *01350000 * POINTS. IN ADDITION THE SLOPE AT THE SECOND SAMPLE *01360000 * POINT OF THE CURRENT GROUP IS FORCED TO BE THE SAME *01370000 * AS THE SLOPE AT THIS POINT CALCULATED USING THE FIT *01380000 * OBTAINED ON THE PREVIOUS GROUP OF POINTS. *01390000 * *01400000 ****************************************************************/01410000 % SKIP(001); 01420000 /* PUT IN CONDITIONS FOR SATISFYING POINT CONDITIONS 01430000 */01440000 % SKIP(001); 01450000 M = 0; 01460000 DO I = INDEX TO INDEX - 1 + NPTS_PER_FIT; 01470000 IF I = INDEX + 1 THEN GO TO PT_CDN_LP; 01480001 M = M + 1; 01490000 X = (I - INDEX)*DX; 01500000 IF M=1 01510000 THEN B(M,1) = F1 - F2 + FP2*DX; 01520000 ELSE B(M,1) = F(I) - F2 - FP2*(X - DX); 01530000 DO J = 1 TO NTERMS - 2; 01540000 A(M,J) = J*DX**(J+1) - (J+1)*DX**J*X + X**(J+1); 01550000 END; 01560000 PT_CDN_LP: 01570001 END; 01580000 % SKIP(001); 01590000 /* NOW PUT IN LEAST SQUARES SLOPE CONDITION. 01600000 */01610000 % SKIP(001); 01620000 M = M + 1; 01630000 B(M,1) = FP1 - FP2; 01640000 DO J = 1 TO NTERMS - 2; 01650000 A(M,J) = -(J+1)*DX**J; 01660000 END; 01670000 CALL MLSQ (A,B,IPTS_PER_FIT,NTERMSM2 ,NUM1); 01680000 % SKIP(001); 01690000 /* RESHUFFLE B TO INCLUDE A1 WHICH HAS BEEN LEFT OUT 01700000 IN ORDER TO FORCE THE SLOPE CONDITION AT X = SECOND 01710000 SAMPLE POINT. 01720000 */01730000 % SKIP(001); 01740000 DUM0 = F2 - FP2*DX; 01750000 DUM = FP2; 01760000 DO J = 1 TO NTERMS - 2; 01770000 B(J+2,1) = B(J,1); 01780000 DUM0 = DUM0 + J*B(J+2,1)*DX**(J+1); 01790000 DUM = DUM-(J+1)* B(J+2,1) * DX **J; 01800000 END; 01810000 B(1,1) = DUM0; 01820000 B(2,1) = DUM; 01830000 % SKIP(001); 01840000 /* STORE NEEDED VALUES FOR NEXT FIT 01850000 */01860000 % SKIP(001); 01870000 F1 = F(INDEX + NPTS_PER_JUMP); 01880000 IF NPTS_PER_JUMP=1 THEN F1 = F2S; 01890000 F2S= F(INDEX + NPTS_PER_JUMP+1); 01900000 XD1 = NPTS_PER_JUMP * DX; 01910000 XD2 = XD1 + DX; 01920000 FP1 = 0.0; 01930000 FP2 = 0.0; 01940000 DO J = 2 TO NTERMS; 01950000 FP1 = FP1 + (J-1)*B(J,1)*XD1**(J-2); 01960000 FP2 = FP2 + (J-1)* B(J,1) * XD2 ** (J-2); 01970000 END; 01980000 END SPLINE_LSTSQRS_FITTER; 01990000 % SKIP(002); 02000000 FREE A, B; 02010000 END LSPYSMO; 02020000 % PAGE; 02030000 MLSQ: 02040000 PROCEDURE (A,B,M,N,K); 02050000 % SKIP(002); 02060000 /*************************************************/ 02070001 /* MLSQ CALCULATES "X" SATISFYING "AX = B", */ 02080001 /* THAT IS, THE SOLUTION OF A SYSTEM OF LINEAR */ 02090001 /* EQUATIONS USING HOUSEHOLDER TRANSFORMATIONS. */ 02100001 /* THE LEAST-SQUARES SOLUTION IS OBTAINED IN */ 02110001 /* CASE OF AN OVERDETERMINED SYSTEM OF EQUATIONS.*/ 02120000 /*************************************************/ 02130001 % SKIP(001); 02140000 DECLARE 02150000 ( A(*,*), B(*,*), PIVR, MAXA) 02160000 BINARY FLOAT, /* SINGLE PRECISION VERSION */ 02170000 (AUX(N),H,SIG,BETA) 02180000 BINARY FLOAT(53), 02190000 (TOL, PIV(N)) 02200000 BINARY FLOAT, 02210000 (I,J,K,L,M,N,PIVI,LM,LN,LK) 02220000 BINARY FIXED; 02230000 LM = M; 02240000 LN = N; 02250000 LK = K; 02260000 SIG = 0; 02270000 IF LM >= LN THEN /* IF M LESS THAN N */ 02280000 IF LN >= 0 THEN /* OR IF N NOT POSITIVE */ 02290000 IF LK > 0 THEN DO; /* OR IF K NOT POSITIVE */ 02300000 /* THEN BYPASS OPERATION */ 02310000 DO L = 1 TO LN; /* CALCULATE SCALARPRODUCTS OF*/ 02320000 H=0; /* COLUMNS. */ 02330000 DO I = 1 TO LM; 02340000 H = H + MULTIPLY(A(I,L),A(I,L),53); 02350000 END; 02360000 IF H >= SIG THEN DO; 02370000 SIG = H; /* SAVE MAX SCALARPRODUCT */ 02380000 PIVI = L; /* SAVE SUBSCRIPT OF PIVOT */ 02390000 END; 02400000 AUX(L),PIV(L) = H; 02410000 END; 02420000 % SKIP(002); 02430000 /* DECOMPOSITION LOOP */ 02440000 % SKIP(001); 02450000 DO L = 1 TO LN; 02460000 TOL = PIV(PIVI); /* ORIGINAL LENGTH OF PIVOT */ 02470000 IF PIVI > L THEN DO; /* SHOULD INTERCHANGE COLUMN */ 02480000 H = AUX(L); /* INTERCHANGE SCALARPRODUCTS */ 02490000 AUX(L) = AUX(PIVI); 02500000 PIV(PIVI) = PIV(L); 02510000 AUX(PIVI) = H; 02520000 DO J = L TO LM; /* INTERCHANGE LOWER PART OF */ 02530000 PIVR = A(J,L); /* COLUMNS OF A */ 02540000 A(J,L) = A(J,PIVI); 02550000 A(J,PIVI) = PIVR; 02560000 END; 02570000 END; 02580000 IF L > 1 THEN DO ; /* RECALCULATE COLUMN LENGTH */ 02590000 SIG = 0; /* TO AVOID ROUND-OFF ERRORS */ 02600000 DO I = L TO LM ; 02610000 SIG = SIG + MULTIPLY(A(I,L),A(I,L),53); 02620000 END; 02630000 END; 02640000 IF TOL = 0 THEN DO; 02650000 TOL = 1; 02660000 END; 02670000 BETA = TOL * 1E-10; 02680000 IF SIG <= BETA THEN DO; /* INDICATE LOSS OF SIGNIFI. */ 02690000 IF SIG <= 0 THEN SIG =BETA; 02700000 END; 02710000 SIG = SQRT(SIG); 02720000 H = A(L,L); 02730000 IF H < 0 THEN SIG = -SIG ; /* FORCE SIGN(SIG) TO SIGN(H)*/ 02740000 PIV(L) = PIVI ; /* SAVE INTERCHANG INFO. */ 02750000 A(L,L),BETA = H + SIG; /* TRANSFORM DIAGONIAL */ 02760000 AUX(L) = -SIG; /* SAVE DIAGONIAL ELEMENT */ 02770000 BETA = SIG * BETA; 02780000 /* TRANSFORM OF SUBMATRIX A */ 02790000 PIVR = 0; 02800000 DO J = L+1 TO LN; 02810000 H =0; 02820000 DO I = L TO LM; 02830000 H = H + MULTIPLY(A(I,L),A(I,J),53); 02840000 END; 02850000 SIG = H/BETA; 02860000 DO I = LM TO L BY -1; 02870000 H = A(I,J); 02880000 A(I,J) = H - A(I,L)*SIG; 02890000 END; 02900000 H = A(L,J); 02910000 H,AUX(J) = AUX(J) - H*H; 02920000 IF H >= PIVR THEN DO; 02930000 PIVR = H; 02940000 PIVI = J; 02950000 END; 02960000 END; 02970000 DO J = 1 TO K; 02980000 H = 0; 02990000 DO I = L TO LM; 03000000 H = H + MULTIPLY(A(I,L),B(I,J),53); 03010000 END; 03020000 MAXA = H/BETA; 03030000 DO I = L TO LM; 03040000 B(I,J) = B(I,J) - A(I,L)*MAXA; 03050000 END; 03060000 END; 03070000 END; 03080000 DO J = LN TO 1 BY -1; 03090000 DO I = 1 TO LK; 03100000 H = B(J,I); 03110000 DO L = J+1 TO LN; 03120000 H = H- MULTIPLY(A(J,L),B(L,I),53); 03130000 END; 03140000 PIVI = PIV(J); 03150000 B(J,I) = B(PIVI,I); 03160000 B(PIVI,I) = H/AUX(J); 03170000 END; 03180000 END; 03190000 IF LN < LM THEN DO J = 1 TO LK; 03200000 H=0; 03210000 DO I = LN + 1 TO LM ; 03220000 H = H + MULTIPLY(B(I,J),B(I,J),53); 03230000 END; 03240000 B(LM,J) = H; 03250000 END; 03260000 END; 03270000 END MLSQ; 03280000