CTITLEMSFIT -- LEAST SQUARES POLYNOMIAL CURVE FIT 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR UNKNOWN 00000020 CA DESIGNER UNKNOWN - PASS ROUTINE CONVERTED BY 00000030 CA GALEN WHIPPLE 08-01-77 00000040 CA LANGUAGE FORTRAN 00000050 CA SYSTEM IBM AND CRAY 00000050 CA WRITTEN MO-DA-YR 00000060 C REVISED MO-DA-YR BY PROGRAMMER FOR REASON. 00000070 C REVISED 00000080 CA 00000090 CA 00000100 CA CALL MSFIT (X, Y, FITORD, COEFF, LAST) 00000110 CA INPUT X = INDEPENDENT VARIABLE IN INPUT FUNCTION. R4 00000120 CA INPUT Y = DEPENDENT VARIABLE IN INPUT FUNCTION. R4 00000130 CA INPUT FITORD = DEGREE OF DESIRED POLYNOMIAL. I4 00000140 CA OUTPUT COEFF = ARRAY OF POLYNOMIAL COEFFICIENTS, R4 00000150 CA CONSTANT TERM IS FIRST. 00000160 CA INPUT LAST = NUMBER OF INPUT POINTS IN (X,Y) FUNCTION. 00000170 CA 00000180 CA 00000185 CA THIS SUBROUTINE COMPUTES THE SET OF POLYNOMIAL COEFFICIENTS FOR 00000190 CA THE POLYNOMIAL THAT BEST FITS THE INPUT DATA IN A LEAST SQUARES 00000200 CA SENSE. THE REQUIRED SET OF LINEAR EQUATIONS IS SET UP AND THEN 00000210 CA SOLVED USING CRAMERS RULE. IF THE MATRIX OF COEFFICIENTS IS 00000220 CA SINGULAR THE COEFF ARRAY WILL BE RETURNED ALL ZEROS. 00000230 CAEND 00000240 C 00000250 SUBROUTINE MSFIT(X,Y,FITORD,COEFF,LAST) 00000260 DOUBLE PRECISION SCRAT 00000270 DIMENSION X(250),Y(250),T1(250),COEFF(10) 00000280 DIMENSION SX(18),SY(10),G(9),SCRAT(100),R(100) 00000290 INTEGER FITORD,T1 00000300 C *** INITIALIZE *** 00000310 DO 10 I=1,9 00000320 COEFF(I) = 0. 00000330 G(I) = 0. 00000340 SX(I) = 0. 00000350 SX(9+I) = 0. 00000360 10 SY(I) = 0. 00000370 COEFF(10) = 0. 00000380 SY(10) = 0. 00000390 M = 0 00000400 C * 00000410 C *** FORM SUMS *** 00000420 C * 00000430 DO 100 I=1,LAST 00000440 C*** IF(T1(I) .LT. 0) GO TO 100 00000450 M = M+1 00000460 G(1) = X(I) 00000470 DO 20 J=2,9 00000480 20 G(J) = G(J-1)*X(I) 00000490 DO 30 J=1,9 00000500 30 SX(J) = SX(J)+G(J) 00000510 DO 40 J=1,9 00000520 40 SX(J+9) = SX(J+9) + G(J)*G(9) 00000530 SY(1) = SY(1)+Y(I) 00000540 DO 60 J=1,FITORD 00000550 60 SY(J+1) = SY(J+1)+G(J)*Y(I) 00000560 100 CONTINUE 00000570 C * 00000580 C *** FORM MATRIX 00000590 C * 00000600 N = FITORD+1 00000610 DO 200 I=1,N 00000620 IF(I .NE. 1) GO TO 110 00000630 R(1) = M 00000640 L = 1 00000650 K = 2 00000660 N1 = FITORD 00000670 110 DO 120 J=L,N1 00000680 R(K) = SX(J) 00000690 120 K = K+1 00000700 IF(I .EQ. 1) GO TO 200 00000710 L = L+1 00000720 200 N1 = N1+1 00000730 C * 00000740 CALL MLSCO (R,COEFF,SY,SCRAT,N) 00000750 RETURN 00000760 END 00000770