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
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
|
C+
C
C SUBROUTINE POLFT1
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 + ...
C
C USAGE
C CALL POLFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR,ARR,IER)
C
C DESCRIPTION OF PARAMETERS
C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y-DATA POINTS
C (OR - IN CASE MODE = 4 - WEIGHTS FOR POINTS)
C NPTS - NO. OF PAIRS OF DATA POINTS
C NTERMS - NO. OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1)
C MODE - DETERMINES METHOD OF WEIGHTING LEAST SQUARES FIT
C 4 (USER DEFINED) WEIGHT(I) = SIGMAY(I)
C 3 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
C 2 (NO WEIGHTING) WEIGHT(I) = 1.
C 1 (STATISTICAL) WEIGHT(I) = 1./Y(I)
C A - ARRAY OF COEFFICIENTS OF POLYNOMIAL
C CHISQR - REDUCED CHISQUARE FOR FIT
C ARR - DOUBLE PRECISION WORK BUFFER; MUST BE AT LEAST
C 400 WORDS LONG IN THE CALLING ROUTINE
C IER - ERROR PARAMETER
C -1 DET=0
C 0 O.K.
C +1 NOT ENOUGH POINTS
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C DET(ARR,NORD) - EVALUATES THE DETERMINANT OF A SYMMETRIC
C TWO DIMENSIONAL MATRIX OF ORDER NORD
C
C COMMENTS
C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C FOR DETAILS OF ALGORITHM SEE "DATA REDUCTION AND ERROR
C ANALYSIS FOR THE PHYSICAL SCIENCES" - PH. R. BEVINGTON
C MC GRAW-HILL BOOK COMPANY
C
C IN THIS SPECIAL VERSION THE ELEMENTS OF COEFFICIENT MATRIX
C ARR ARE NORMALIZED WITH RESPECT TO A VALUE COMPUTED VIA
C LOGARITHMIC INTERPOLATION BETWEEN THE SMALLEST AND LARGEST
C MATRIX ELEMENT
C
C MODIFICATIONS BY A. SCHERMANN - INST. FOR ASTRONOMY,
C UNIV. OF VIENNA, AUSTRIA
C
C Modified to assume x-variable is index of Y array
C-
SUBROUTINE POLFT1(Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR,ARR,IER)
DOUBLE PRECISION SUMX(19),SUMY(10),ARR(10,10),XTERM,YTERM,
1 CHISQ
DIMENSION Y(1),SIGMAY(1),A(1)
IER=0
C
C CHECK DEGREE OF FREEDOM
FREE=NPTS-NTERMS
IF(FREE.GT.0.) GOTO 11
IER=1
GOTO 80
C
C ACCUMULATE WEIGHTED SUMS
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=I
YI=Y(I)
31 GOTO(32,37,38,39),MODE
32 IF(YI)35,37,33
33 WEIGHT=1./YI
GOTO 41
35 WEIGHT=-1./YI
GOTO 41
37 WEIGHT=1.
GOTO 41
38 WEIGHT=1./SIGMAY(I)**2
GOTO 41
39 WEIGHT=SIGMAY(I)
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 GET LARGEST AND SMALLEST MATRIX ELEMENT (FOR NORMALIZATION)
XTERM=SUMX(1)
YTERM=XTERM
DO 100 I=2,NMAX
IF(SUMX(I).GT.XTERM) XTERM=SUMX(I)
IF(SUMX(I).LT.YTERM) YTERM=SUMX(I)
100 CONTINUE
DO 110 I=1,NTERMS
IF(SUMY(I).GT.XTERM) XTERM=SUMY(I)
IF(SUMY(I).LT.YTERM) YTERM=SUMY(I)
110 CONTINUE
IF(YTERM.LE.0.) YTERM=1.D0
C
C LOGARITHMIC INTERPOLATION OF NORMALIZATION VALUE
XTERM=1.D1**((DLOG10(XTERM)+DLOG10(YTERM))/2.)
C
C CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS
51 DO 54 J=1,NTERMS
DO 54 K=1,NTERMS
N=J+K-1
54 ARR(J,K)=SUMX(N)/XTERM
DELTA=DET(ARR,NTERMS)
IF(DELTA.NE.0.) GOTO 61
57 CHISQR=0.
DO 59 J=1,NTERMS
59 A(J)=0.
IER=-1
GOTO 80
61 DO 70 L=1,NTERMS
62 DO 66 J=1,NTERMS
DO 65 K=1,NTERMS
N=J+K-1
65 ARR(J,K)=SUMX(N)/XTERM
66 ARR(J,L)=SUMY(J)/XTERM
70 A(L)=DET(ARR,NTERMS)/DELTA
C
C CALCULATE CHISQUARE
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)
77 CHISQR=CHISQ/FREE
80 RETURN
END
c-----------------------------------------------------------
FUNCTION POLYNO (COE,NPOL,IX)
C
C EVALUATE A POLYNOMIAL OF ORDER NPOL WITH COEFFICIENTS COE(1),...,
C COE (NPOL+1) FOR AN INDEX IX
C
DIMENSION COE(*)
C
IF(NPOL.GT.0) GO TO 10
POLYNO=COE(1)
RETURN
C
10 POLYNO=COE(NPOL+1)
DO 20 I=NPOL,1,-1
20 POLYNO=POLYNO*IX+COE(I)
C
RETURN
END
c-----------------------------------------------------------
C+
C
C FUNCTION DET
C
C PURPOSE
C CALCULATE DETERMINANT OF A SQUARE MATRIX
C
C USAGE
C DET = DET (ARR,NORDER)
C
C DESCRIPTION OF PARAMETERS
C ARRAY - MATRIX
C NORDER - DEGREE OF MATRIX
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C COMMENTS
C THIS SUBPROGRAM DESTROYS THE INPUT MATRIX ARRAY
C
FUNCTION DET(ARRAY,NORDER)
DOUBLE PRECISION ARRAY(10,10),SAVE
DET=1.
DO 90 K=1,NORDER
C INTERCHANGE COLUMNS, IF DIAGONAL ELEMENT IS ZERO
IF(ARRAY(K,K).NE.0.) GOTO 40
DO 10 J=K,NORDER
IF(ARRAY(K,J).NE.0.) GOTO 20
10 CONTINUE
DET=0.
GOTO 100
20 DO 30 I=K,NORDER
SAVE=ARRAY(I,J)
ARRAY(I,J)=ARRAY(I,K)
30 ARRAY(I,K)=SAVE
DET=-DET
C SUBTRACT ROW K FROM LOWER ROWS TO GET DIAGONAL MATRIX
40 DET=DET*ARRAY(K,K)
IF(K-NORDER.GE.0) GOTO 90
K1=K+1
DO 50 I=K1,NORDER
DO 50 J=K1,NORDER
50 ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
90 CONTINUE
100 RETURN
END
|
