CTITLEUSPFIT -- UTILITY POLYNOMIAL FIT ROUTINE 00010012 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR RANDY PROBST 00020012 CA DESIGNER RANDY PROBST 00030012 CA LANGUAGE FORTRAN 00040012 CA SYSTEM IBM AND CRAY 00050012 CA WRITTEN 09-19-86 00060012 C REVISED 00070012 CA 00080012 CA CALL USPFIT (X,Y,INC,NSMP,LCPI,PASS,ERR,B,PRTFLG,IPR) 00090012 CA IN X = ARRAY OF X-COORDINATES R4 00100012 CA IN Y = ARRAY OF Y-COORDINATES R4 00110012 CA IN INC = INDEX INCREMENT FOR X&Y ARRAYS I4 00120012 CA IN NSMP = NUMBER OF ELEMENTS OF X&Y ARRAYS I4 00130012 CA IN LCPI = SAMPLE INTERVAL I4 00140012 CA IN PASS = NUMBER OF PASSES TO MAKE ON FIT I4 00150012 CA OUT ERR = RMS ERROR ON FIT R4 00160012 CA OUT B = ARRAY CONTAINING POLYNOMIAL CONSTANTS R4 00170012 CA A,B,C => (Y = A + BX + CXX) 00180012 CA IN PRTFLG = DIAGNOSTIC PRINT FLAG (=64 CAUSES PRINT) I4 00190012 CA IN IPR = LOGICAL UNIT NUMBER FOR PRINT I4 00200012 CA 00210012 CA 00220012 CA THIS ROUTINE TAKES AN ARRAY OF X COORDINATES AND Y COORDINATES 00230012 CA AND PERFORMS A POLYNOMIAL FIT CALCULATING THE CONSTANTS IN THE 00240012 CA EXPRESSION (Y = A + BX +CXX). THE VALUES OF A, B, & C ARE 00250012 CA RETURNED IN THE ARRAY B. MULTIPLE PASSES MAY BE MADE THROWING 00260012 CA OUT VALUES OF MORE THAN TWO RMS VARIATIONS FROM THE CURVE WITH 00270012 CA 500 MILLISECONDS BEING A MAXIMUM. 00280012 CA 00290012 SUBROUTINE USPFIT (X,Y,INC,NSMP,LCPI,PASS,ERR,B,PRTFLG,IPR) 00300008 IMPLICIT INTEGER (A-Z) 00310001 C 00320001 INTEGER X(1), Y(1) 00330001 C 00340001 REAL X1, X2, X3, X4, XY, X2Y, Y1, TMP, THR 00350005 REAL ERR 00360001 C 00370001 REAL A(3,3), B(3) 00380001 C 00390001 LOOP = 1 00400001 ERR = 99999.0 00410001 B(1) = 0.0 00420001 B(2) = 0.0 00430001 B(3) = 0.0 00440001 C 00450001 50 CONTINUE 00460001 X1 = 0.0 00470001 Y1 = 0.0 00480001 XY = 0.0 00490001 X2 = 0.0 00500001 X2Y = 0.0 00510001 X3 = 0.0 00520001 X4 = 0.0 00530001 M = 0 00540001 K = 0 00550001 C 00560001 DO 100 J = 1, NSMP 00570001 I = (J-1) * INC + 1 00580001 IF (Y(I) .LE. 0) GO TO 100 00590001 IF (LOOP .GT. 1) THEN 00600001 TMP = B(1) + B(2)*X(I) + B(3)*X(I)*X(I) 00610001 IF (ABS(TMP-Y(I)) .GT. THR) THEN 00620001 K = K + 1 00630001 GO TO 100 00640001 END IF 00650001 END IF 00660001 X1 = X1 + X(I) 00670001 Y1 = Y1 + Y(I) 00680001 XY = XY + X(I) * Y(I) 00690001 TMP = X(I) * X(I) 00700001 X2 = X2 + TMP 00710001 X2Y = X2Y + TMP * Y(I) 00720001 TMP = TMP * X(I) 00730001 X3 = X3 + TMP 00740001 X4 = X4 + TMP * X(I) 00750001 M = M + 1 00760001 100 CONTINUE 00770001 C 00780001 C CHECK FOR NO IMPROVEMENTS TO ANOTHER PASS 00790001 IF (K .EQ. 0 .AND. LOOP .GT. 1) GO TO 400 00800001 C 00810001 C CHECK FOR NO VALUES 00820001 IF (M .LT. 3) GO TO 400 00830014 C 00840001 A(1,1) = M 00850001 A(2,1) = X1 00860001 A(3,1) = X2 00870001 A(1,2) = X1 00880001 A(2,2) = X2 00890001 A(3,2) = X3 00900001 A(1,3) = X2 00910001 A(2,3) = X3 00920001 A(3,3) = X4 00930001 B(1) = Y1 00940001 B(2) = XY 00950001 B(3) = X2Y 00960001 C 00970001 CALL MATINV (A, 3, 3, B, 1) 00980001 ERR = 0.0 00990001 DO 200 J = 1, NSMP 01000001 I = (J-1) * INC + 1 01010001 IF (Y(I) .LE. 0) GO TO 200 01020001 TMP = B(1) + B(2)*X(I) + B(3)*X(I)*X(I) 01030001 IF (LOOP .GT. 1 .AND. ABS(TMP-Y(I)) .GT. THR) GO TO 200 01040001 ERR = ERR + (TMP - Y(I)) * (TMP - Y(I)) 01050001 200 CONTINUE 01060001 ERR = SQRT(ERR/M) 01070001 THR = ERR * LCPI * 2.0 01080010 IF (THR .GT. 500.0) THR = 500.0 01090011 THR = THR / LCPI 01100010 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 01110001 IF(IAND(PRTFLG,64).NE.0) 01120001 + WRITE(IPR,2000) (B(I),I=1,3),ERR*LCPI,THR*LCPI,LOOP,NSMP,M 01130010 2000 FORMAT(' USPFIT -A,B,C,ERR,THR,LOOP,NSMP,M-', 01140008 + 3E13.4, E21.10, F8.2, 3I8) 01150008 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 01160001 C 01170001 LOOP = LOOP + 1 01180001 IF (PASS .GE. LOOP) GO TO 50 01190001 C 01200001 400 CONTINUE 01210001 C 01220001 RETURN 01230001 END 01240001