Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap1_198111 - decus/20-0020/lscf.sta
There are 2 other files named lscf.sta in the archive. Click here to see a list.
5000'  NAME--LSCF
5010'
5020'  DESCRIPTION--LEAST SQUARES POLYNOMIAL CURVE FIT SUBROUTINE
5030'  FITS A SET OF X,Y PAIRS TO THE EQUATION
5040'      Y = A + B*X + C*X 2 + D*X 3 + .......
5050'
5060'  SOURCE--L.C.SEMPREBON, RADIOPHYSICS LABORATORY,
5070'  DARTMOUTH COLLEGE, HANOVER, N.H. 03755
5080'
5090'  INSTRUCTIONS--BEFORE ENTERING THE SUBROUTINE THE FOLLOWING
5100'  MUST BE DEFINED:
5110'
5120'  N9      = NUMBER OF X,Y PAIRS.
5130'  D(N9,3) = THE DATA ARRAY. COLUMN 1 CONTAINS X VALUES,
5140'            COLUMS 2 CONTAINS Y VALUES, COLUMN 3 WILL CONTAIN
5150'            THE FITTED Y VALUES---ONLY COLUMNS 1 AND 2 NEED TO
5160'            BE DEFINED.
5170'  D9      = DESIRED DEGREE OF THE POLYNOMIAL.
5180'
5190'  DIM STATEMENTS FOR THE FOLLOWING MATRICES MUST BE GIVEN, WITH
5200'  "*" REPLACED WITH THE MAXIMUM VALUE EXPECTED FOR D9+1
5210'  DIM V(*,*),W(*),X(*,*),Y(*),Z(2*),D(N9,3)
5220'
5230'  ALL VARIABLES USED IN THE SUBROUTINE END IN 9.
5240'
5250'  OUTPUT:Y(D9)  = SOLUTION VECTOR, WHERE A=Y(0),B=Y(1), ETC.
5260'           S9     = SUM OF THE SQUARES OF THE DIFFERENCES
5270'                    BETWEEN THE COMPUTED AND INITIAL VALUES OF Y.
5280'           D(N9,2) = COLUMN 2 OF THE INITIALIZED DATA ARRAY
5290'                    WILL CONTAIN THE FITTED VALUES OF Y.
5300'
5310'
5320'  *  *  *  *  *  *   SUBROUTINE   *  *  *  *  *  *  *  *  *
5330'
5340 MAT Z=ZER(2*D9+1)
5350 MAT W=ZER(D9+1)
5360 LET Z(1)=N9
5370 FOR I9=1 TO N9
5380 LET X9=D(I9,1)
5390 LET Y9=D(I9,2)
5400 FOR J9=2 TO 2*D9+1
5410 LET Z(J9)=Z(J9)+X9
5420 LET X9=X9*D(I9,1)
5430 NEXT J9
5440 LET X9=1
5450 FOR J9=1 TO D9+1
5460 LET W(J9)=W(J9)+X9*Y9
5470 LET X9=X9*D(I9,1)
5480 NEXT J9
5490 NEXT I9
5500 MAT V=ZER(D9+1,D9+1)
5510 FOR I9=1 TO D9+1
5520 FOR J9=1 TO D9+1
5530 LET V(I9,J9)=Z(I9+J9-1)
5540 NEXT J9
5550 NEXT I9
5560 MAT X=INV(V)
5570 MAT Y=X*W
5580 LET S9=0
5590 FOR I9=1 TO N9
5600 LET X9=1
5610 LET D(I9,3)=0
5620 FOR J9=1 TO D9+1
5630 LET D(I9,3)=D(I9,3)+X9*Y(J9)
5640 LET X9=D(I9,1)*X9
5650 NEXT J9
5660 LET S9=S9+(D(I9,2)-D(I9,3))^2
5670 NEXT I9
5680 RETURN