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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
c subroutine legfit.f
c
c source
c Bevington, pages 155-157.
c
c purpose
c make a least-squares fit to data with a legendre polynomial
c y = a(1) + a(2)*x + a(3)*(3x**2-1)/2 + ...
c = a(1) * (1. + b(2)*x + b(3)*(3x**2-1)/2 + ... )
c where x = cos(theta)
c
c usage
c call legfit (theta, y, sigmay, npts, norder, neven, mode, ftest,
c yfit, a, sigmaa, b, sigmab, chisqr)
c
c description of parameters
c theta - array of angles (in degrees) of the data points
c y - array of data points for dependent variable
c sigmay - array of standard deivations for y data points
c npts - number of pairs of data points
c norder - highest order of polynomial (number of terms - 1)
c neven - determines odd or even character of polynomial
c +1 fits only to even terms
c 0 fits to all terms
c -1 fits only to odd terms (plus constant term)
c mode - determines mode 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 ftest - array of values of f(l) for an f test
c yfit - array of calculated values of y
c a - array of coefficients of polynomial
c sigmaa - array of standard deviations for coefficients
c b - array of normalized relative coefficients
c sigmab - array of standard deviations for relative coefficients
c chisqr - reduced chi square for fit
c
c subroutines and function subprograms required
c matinv (array, nterms, det)
c inverts a symmetric two-dimensional matrix of degree nterms
c and calculates its determinant
c
c comments
c dimension statement valid for npts up to 100 and order up to 9
c dcos changed to cos in statement 31
c
subroutine legfit (theta,y,sigmay,npts,norder,neven,mode,
*ftest,yfit,a,sigmaa,b,sigmab,chisqr)
double precision cosine,p,beta,alpha,chisq
dimension theta(1),y(1),sigmay(1),ftest(1),yfit(1),
*a(1),sigmaa(1),b(1),sigmab(1)
dimension weight(100),p(100,10),beta(10),alpha(10,10)
c
c accumulate weights and legendre polynomials
c
11 nterms=1
ncoeff=1
jmax=norder+1
20 do 40 i=1,npts
21 if (mode) 22,27,29
22 if (y(i)) 25,27,23
23 weight(i)=1./y(i)
goto 31
25 weight(i)=1./(-y(i))
goto 31
27 weight(i)=1.
goto 31
29 weight(i)=1./sigmay(i)**2
31 cosine=cos(0.01745329252*theta(i))
p(i,1)=1.
p(i,2)=cosine
do 36 l=2,norder
fl=l
36 p(i,l+1)=((2.*fl-1.)*cosine*p(i,l)-(fl-1.)*p(i,l-1))/fl
40 continue
c
c accumulate matrices alpha and beta
c
51 do 54 j=1,nterms
beta(j)=0.
do 54 k=1,nterms
54 alpha(j,k)=0.
61 do 66 i=1,npts
do 66 j=1,nterms
beta(j)=beta(j)+p(i,j)*y(i)*weight(i)
do 66 k=j,nterms
alpha(j,k)=alpha(j,k)+p(i,j)*p(i,k)*weight(i)
66 alpha(k,j)=alpha(j,k)
c
c delete fixed coefficients
c
70 if (neven) 71,91,81
71 do 76 j=3,nterms,2
beta(j)=0.
do 75 k=1,nterms
alpha(j,k)=0.
75 alpha(k,j)=0.
76 alpha(j,j)=1.
goto 91
81 do 86 j=2,nterms,2
beta(j)=0.
do 85 k=1,nterms
alpha(j,k)=0.
85 alpha(k,j)=0.
86 alpha(j,j)=1.
c
c invert curvature matrix alpha
c
91 do 95 j=1,jmax
a(j)=0.
sigmaa(j)=0.
b(j)=0.
95 sigmab(j)=0.
do 97 i=1,npts
97 yfit(i)=0.
101 call matinv (alpha,nterms,det)
if (det) 111,103,111
103 chisqr=0.
goto 170
c
c calculate coefficients, fit, and chi square
c
111 do 115 j=1,nterms
do 113 k=1,nterms
113 a(j)=a(j)+beta(k)*alpha(j,k)
do 115 i=1,npts
115 yfit(i)=yfit(i)+a(j)*p(i,j)
121 chisq=0.
do 123 i=1,npts
123 chisq=chisq+(y(i)-yfit(i))**2*weight(i)
free=npts-ncoeff
chisqr=chisq/free
c
c test for end of fit
c
131 if (nterms-jmax) 132,151,151
132 if (ncoeff-2) 133,134,141
133 if (neven) 137,137,135
134 if (neven) 135,137,135
135 nterms=nterms+2
goto 138
137 nterms=nterms+1
138 ncoeff=ncoeff+1
chisq1=chisq
goto 51
141 fvalue=(chisq1-chisq)/chisqr
if (ftest(nterms)-fvalue) 134,143,143
143 if (neven) 144,146,144
144 nterms=nterms-2
goto 147
146 nterms=nterms-1
147 ncoeff=ncoeff-1
jmax=nterms
149 goto 51
c
c calculate remainder of output
c
151 if (mode) 152,154,152
152 varnce=1.
goto 155
154 varnce=chisqr
155 do 156 j=1,nterms
156 sigmaa(j)=dsqrt(varnce*alpha(j,j))
161 if (a(1)) 162,170,162
162 do 166 j=2,nterms
if (a(j)) 164,166,164
164 b(j)=a(j)/a(1)
165 sigmab(j)=b(j)*dsqrt((sigmaa(j)/a(j))**2+(sigmaa(1)/a(1))**2
*-2.*varnce*alpha(j,1)/(a(j)*a(1)))
166 continue
b(1)=1.
170 return
end
|
