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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/curfit.h>
include "curfitdef.h"
# CVACCUM -- Procedure to add a data point to the set of normal equations.
# The inner products of the basis functions are added into the CV_ORDER(cv)
# by CV_NCOEFF(cv) array MATRIX. The first row of MATRIX
# contains the main diagonal of the matrix followed by
# the CV_ORDER(cv) lower diagonals. This method of storing MATRIX
# minimizes the space required by large symmetric, banded matrices.
# The inner products of the basis functions and the data ordinates are
# stored in VECTOR which has CV_NCOEFF(cv) elements. The integers left
# and leftm1 are used to determine which elements of MATRIX and VECTOR
# are to receive the data.
procedure cvaccum (cv, x, y, w, wtflag)
pointer cv # curve descriptor
real x # x value
real y # y value
real w # weight of the data point
int wtflag # type of weighting desired
int left, i, ii, j
real bw
pointer xzptr
pointer mzptr, mzzptr
pointer vzptr
begin
# increment number of points
CV_NPTS(cv) = CV_NPTS(cv) + 1
# calculate the weights
switch (wtflag) {
case WTS_UNIFORM, WTS_SPACING:
w = 1.0
case WTS_USER:
# user defined weights
case WTS_CHISQ:
# data assumed to be scaled to photons with Poisson statistics
if (y > 0.0)
w = 1.0 / y
else if (y < 0.0)
w = - 1.0 / y
else
w = 0.0
default:
w = 1.0
}
# calculate all non-zero basis functions for a given data point
switch (CV_TYPE(cv)) {
case CHEBYSHEV:
left = 0
call rcv_b1cheb (x, CV_ORDER(cv), CV_MAXMIN(cv), CV_RANGE(cv),
XBASIS(CV_XBASIS(cv)))
case LEGENDRE:
left = 0
call rcv_b1leg (x, CV_ORDER(cv), CV_MAXMIN(cv), CV_RANGE(cv),
XBASIS(CV_XBASIS(cv)))
case SPLINE3:
call rcv_b1spline3 (x, CV_NPIECES(cv), -CV_XMIN(cv),
CV_SPACING(cv), XBASIS(CV_XBASIS(cv)), left)
case SPLINE1:
call rcv_b1spline1 (x, CV_NPIECES(cv), -CV_XMIN(cv),
CV_SPACING(cv), XBASIS(CV_XBASIS(cv)), left)
case USERFNC:
left = 0
call rcv_b1user (cv, x)
}
# index the pointers
xzptr = CV_XBASIS(cv) - 1
vzptr = CV_VECTOR(cv) + left - 1
mzptr = CV_MATRIX(cv) + CV_ORDER(CV) * (left - 1)
# accumulate the data point into the matrix and vector
do i = 1, CV_ORDER(cv) {
# calculate the non-zero basis functions
bw = XBASIS(xzptr+i) * w
# add the inner product of the basis functions and the ordinate
# into the appropriate element of VECTOR
VECTOR(vzptr+i) = VECTOR(vzptr+i) + bw * y
# accumulate the inner products of the basis functions into
# the apprpriate element of MATRIX
mzzptr = mzptr + i * CV_ORDER(cv)
ii = 0
do j = i, CV_ORDER(cv) {
MATRIX(mzzptr+ii) = MATRIX(mzzptr+ii) + XBASIS(xzptr+j) * bw
ii = ii + 1
}
}
end
|
