1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
C SUBROUTINE POLFIT.F
C
C SOURCE
C BEVINGTON, PAGES 140-142.
C
C PURPOSE
C MAKE A LEAST-SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE
C Y = A(1) + A(2)*X + A(3)*X**2 + A(3)*X**3 + . . .
C
C USAGE
C CALL POLFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, CHISQR)
C
C DESCRIPTION OF PARAMETERS
C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS
C NPTS - NUMBER OF PAIRS OF DATA POINTS
C NTERMS - NUMBER OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1)
C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
C +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
C 0 (NO WEIGHTING) WEIGHT(I) = 1.
C -1 (STATISTICAL) WEIGHT(I) = 1./Y(I)
C A - ARRAY OF COEFFICIENTS OF POLYNOMIAL
C CHISQR - REDUCED CHI SQUARE FOR FIT
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C DETERM (ARRAY, NORDER)
C EVALUATES THE DETERMINANT OF A SYMETRICAL
C TWO-DIMENSIONAL MATRIX OF ORDER NORDER
C
C COMMENTS
C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C
SUBROUTINE POLFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR)
DOUBLE PRECISION SUMX,SUMY,XTERM,YTERM,ARRAY,CHISQ
DIMENSION X(1),Y(1),SIGMAY(1),A(1)
DIMENSION SUMX(19),SUMY(10),ARRAY(10,10)
C
C ACCUMULATE WEIGHTED SUMS
C
11 NMAX=2*NTERMS-1
DO 13 N=1,NMAX
13 SUMX(N)=0.
DO 15 J=1,NTERMS
15 SUMY(J)=0.
CHISQ=0.
21 DO 50 I=1,NPTS
XI=X(I)
YI=Y(I)
31 IF (MODE) 32,37,39
32 IF (YI) 35,37,33
33 WEIGHT=1./YI
GOTO 41
35 WEIGHT=1./(-YI)
GOTO 41
37 WEIGHT=1.
GOTO 41
39 WEIGHT=1./SIGMAY(I)**2
41 XTERM=WEIGHT
DO 44 N=1,NMAX
SUMX(N)=SUMX(N)+XTERM
44 XTERM=XTERM*XI
45 YTERM=WEIGHT*YI
DO 48 N=1,NTERMS
SUMY(N)=SUMY(N)+YTERM
48 YTERM=YTERM*XI
49 CHISQ=CHISQ+WEIGHT*YI**2
50 CONTINUE
C
C CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS
C
51 DO 54 J=1,NTERMS
DO 54 K=1,NTERMS
N=J+K-1
54 ARRAY(J,K)=SUMX(N)
DELTA=DETERM (ARRAY,NTERMS)
IF (DELTA) 61,57,61
57 CHISQR=0.
DO 59 J=1,NTERMS
59 A(J)=0.
GOTO 80
61 DO 70 L=1,NTERMS
62 DO 66 J=1,NTERMS
DO 65 K=1,NTERMS
N=J+K-1
65 ARRAY(J,K)=SUMX(N)
66 ARRAY(J,L)=SUMY(J)
70 A(L)=DETERM (ARRAY,NTERMS) /DELTA
C
C CALCULATE CHI SQUARE
C
71 DO 75 J=1,NTERMS
CHISQ=CHISQ-2.*A(J)*SUMY(J)
DO 75 K=1,NTERMS
N=J+K-1
75 CHISQ=CHISQ+A(J)*A(K)*SUMX(N)
76 FREE=NPTS-NTERMS
77 CHISQR=CHISQ/FREE
80 RETURN
END
|
