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
|
C SUBROUTINE CHIFIT.F
C
C SOURCE
C BEVINGTON, PAGES 228-231.
C
C PURPOSE
C MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION
C WITH A PARABOLIC EXPANSION OF CHI SQUARE
C
C USAGE
C CALL CHIFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA,
C SIGMAA, YFIT, 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 PARAMETERS
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 PARAMETERS
C DELTAA - ARRAY OF INCREMENTS FOR PARAMETERS A
C SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A
C YFIT - ARRAY OF CALCULATED VALUES OF Y
C CHISQR - REDUCED CHI SQUARE FOR FIT
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C FUNCTN (X, I, A)
C EVALUATES THE FITTING FUNCTION FOR THE ITH TERM
C FCHISQ (Y, SIGMAY, NPTS, NFREE, MODE, YFIT)
C EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA
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 NTERMS UP TO 10
C
SUBROUTINE CHIFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA,
*SIGMAA,YFIT,CHISQR)
DOUBLE PRECISION ALPHA
DIMENSION X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),SIGMAA(1),
*YFIT(1)
DIMENSION ALPHA(10,10),BETA(10),DA(10)
REAL FUNCTN
EXTERNAL FUNCTN
C
11 NFREE=NPTS-NTERMS
FREE=NFREE
IF (NFREE) 14,14,16
14 CHISQR=0.
GOTO 120
16 DO 17 I=1,NPTS
17 YFIT(I)=FUNCTN (X,I,A)
CHISQ1=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
C
C EVALUATE ALPHA AND BETA MATRICES
C
20 DO 60 J=1,NTERMS
C
C A(J) + DELTAA(J)
C
21 AJ=A(J)
A(J)=AJ+DELTAA(J)
DO 24 I=1,NPTS
24 YFIT(I)=FUNCTN (X,I,A)
CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
ALPHA(J,J)=CHISQ2-2.*CHISQ1
BETA(J)=-CHISQ2
31 DO 50 K=1,NTERMS
IF (K-J) 33,50,36
33 ALPHA(K,J)=(ALPHA(K,J)-CHISQ2)/2.
ALPHA(J,K)=ALPHA(K,J)
GOTO 50
36 ALPHA(J,K)=CHISQ1-CHISQ2
C
C A(J) + DELTAA(J) AND A(K) + DELTAA(K)
C
41 AK=A(K)
A(K)=AK+DELTAA(K)
DO 44 I=1,NPTS
44 YFIT(I)=FUNCTN (X,I,A)
CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
ALPHA(J,K)=ALPHA(J,K)+CHISQ3
A(K)=AK
50 CONTINUE
C
C A(J) - DELTAA(J)
C
51 A(J)=AJ-DELTAA(J)
DO 53 I=1,NPTS
53 YFIT(I)=FUNCTN (X,I,A)
CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
A(J)=AJ
ALPHA(J,J)=(ALPHA(J,J)+CHISQ3)/2.
BETA(J)=(BETA(J)+CHISQ3)/4.
60 CONTINUE
C
C ELIMINATE NEGATIVE CURVATURE
C
61 DO 70 J=1,NTERMS
IF (ALPHA(J,J)) 63,65,70
63 ALPHA(J,J)=-ALPHA(J,J)
GOTO 66
65 ALPHA(J,J)=0.01
C Changed from DO 70
66 DO 700 K=1,NTERMS
C Changed from 68, 70, 68
IF (K-J) 68,700,68
68 ALPHA(J,K)=0.
ALPHA(K,J)=0.
C New continuation statement
700 CONTINUE
70 CONTINUE
C
C INVERT MATRIX AND EVALUATE PARAMETER INCREMENTS
C
71 CALL MATINV (ALPHA,NTERMS,DET)
DO 76 J=1,NTERMS
DA(J)=0.
74 DO 75 K=1,NTERMS
75 DA(J)=DA(J)+BETA(K)*ALPHA(J,K)
76 DA(J)=0.2*DA(J)*DELTAA(J)
C
C MAKE SURE CHI SQUARE DECREASES
C
81 DO 82 J=1,NTERMS
82 A(J)=A(J)+DA(J)
83 DO 84 I=1,NPTS
84 YFIT(I)=FUNCTN (X,I,A)
CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
IF (CHISQ1-CHISQ2) 87,91,91
87 DO 89 J=1,NTERMS
DA(J)=DA(J)/2.
89 A(J)=A(J)-DA(J)
GOTO 83
C
C INCREMENT PARAMETERS UNTIL CHI SQUARE STARTS TO INCREASE
C
91 DO 92 J=1,NTERMS
92 A(J)=A(J)+DA(J)
DO 94 I=1,NPTS
94 YFIT(I)=FUNCTN (X,I,A)
CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
IF (CHISQ3-CHISQ2) 97,101,101
97 CHISQ1=CHISQ2
CHISQ2=CHISQ3
99 GOTO 91
C
C FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS
C
101 DELTA=1./(1.+(CHISQ1-CHISQ2)/(CHISQ3-CHISQ2))+0.5
DO 104 J=1,NTERMS
A(J)=A(J)-DELTA*DA(J)
104 SIGMAA(J)=DELTAA(J)*SQRT(FREE*ALPHA(J,J))
DO 106 I=1,NPTS
106 YFIT(I)=FUNCTN (X,I,A)
CHISQR=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
111 IF (CHISQ2-CHISQR) 112,120,120
112 DO 113 J=1,NTERMS
113 A(J)=A(J)+(DELTA-1.)*DA(J)
DO 115 I=1,NPTS
115 YFIT(I)=FUNCTN (X,I,A)
CHISQR=CHISQ2
120 RETURN
END
|
