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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/surfit.h>
include "surfitdef.h"
# ISSOLVE -- Procedure to solve the x coefficient matrix for the surface
# coefficients. The inner products of the y basis functions are accumulated
# and stored in the SF_YORDER(sf) by SF_NYCOEFF(sf) array YMATRIX. The
# main diagonal of YMATRIX is stored in the first row of YMATRIX followed
# by the remaining non-zero lower diagonals. The Cholesky factorization
# of YMATRIX is calculated and stored on top of YMATRIX destroying the
# original data. The inner products of the y basis functions and
# and the i-th column of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF
# containing the i-th x coefficients for each line are calculated and
# placed in the i-th row of the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix
# COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF is solved to determine
# the SF_NXCOEFF(sf) by SF_NYCOEFF(sf) surface coefficients. After a
# call to SIFSOLVE the coefficient of the i-th x basis function and the
# j-th y basis function will be found in the j-th column and i-th row
# of COEFF.
procedure issolve (sf, lines, nlines, ier)
pointer sf # pointer to the curve descriptor structure
int lines[ARB] # line numbers included in the fit
int nlines # number of lines fit
int ier # error code
int i, ii, j, k, nxcoeff
pointer ybzptr, ybptr
pointer ylzptr
pointer ymzptr, ymindex
pointer xczptr, xcptr, xcindex
pointer czptr, cptr
pointer left, tleft, rows
pointer sp
begin
# define pointers
ybzptr = SF_YBASIS(sf) - 1
ymzptr = SF_YMATRIX(sf)
xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1
czptr = SF_COEFF(sf) - 1
# zero out coefficient matrix and the y coefficient matrix
call aclrr (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf) * SF_NYCOEFF(sf))
call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf))
# increment the number of points
SF_NYPTS(sf) = nlines
switch (SF_TYPE(sf)) {
case SF_LEGENDRE, SF_CHEBYSHEV:
# accumulate the y value in the y matrix
nxcoeff = SF_NXCOEFF(sf)
do i = 1, SF_YORDER(sf) {
cptr = czptr + i
do k = 1, nxcoeff {
xcptr = xczptr + k
do j = 1, nlines {
xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
COEFF(cptr) = COEFF(cptr) +
YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex)
}
cptr = cptr + SF_NYCOEFF(sf)
}
ii = 0
ybptr = ybzptr
do k = i, SF_YORDER(sf) {
ymindex = ymzptr + ii
do j = 1, nlines
YMATRIX(ymindex) = YMATRIX(ymindex) +
YBASIS(ybzptr+lines[j]) *
YBASIS(ybptr+lines[j])
ii = ii + 1
ybptr = ybptr + SF_NLINES(sf)
}
if (SF_XTERMS(sf) == NO)
nxcoeff = 1
ybzptr = ybzptr + SF_NLINES(sf)
ymzptr = ymzptr + SF_YORDER(sf)
}
case SF_SPLINE3, SF_SPLINE1:
call smark (sp)
call salloc (left, nlines, TY_INT)
call salloc (tleft, nlines, TY_INT)
call salloc (rows, nlines, TY_INT)
ylzptr = SF_YLEFT(sf) - 1
do j = 1, nlines
Memi[left+j-1] = YLEFT(ylzptr+lines[j])
call amulki (Memi[left], SF_YORDER(sf), Memi[rows], nlines)
call aaddki (Memi[rows], SF_YMATRIX(sf), Memi[rows], nlines)
call aaddki (Memi[left], czptr, Memi[left], nlines)
# accumulate the y value in the y matrix
nxcoeff = SF_NXCOEFF(sf)
do i = 1, SF_YORDER(sf) {
call aaddki (Memi[left], i, Memi[tleft], nlines)
do k = 1, nxcoeff {
xcptr = xczptr + k
do j = 1, nlines {
cptr = Memi[tleft+j-1]
xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
COEFF(cptr) = COEFF(cptr) + YBASIS(ybzptr+lines[j]) *
XCOEFF(xcindex)
}
call aaddki (Memi[tleft], SF_NYCOEFF(sf), Memi[tleft],
nlines)
}
ii = 0
ybptr = ybzptr
do k = i, SF_YORDER(sf) {
do j = 1, nlines {
ymindex = Memi[rows+j-1] + ii
YMATRIX(ymindex) = YMATRIX(ymindex) +
YBASIS(ybzptr+lines[j]) *
YBASIS(ybptr+lines[j])
}
ii = ii + 1
ybptr = ybptr + SF_NLINES(sf)
}
ybzptr = ybzptr + SF_NLINES(sf)
call aaddki (Memi[rows], SF_YORDER(sf), Memi[rows], nlines)
}
call sfree (sp)
}
# return if not enough points
ier = OK
if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) {
ier = NO_DEG_FREEDOM
return
}
# calculate the Cholesky factorization of the y matrix
call sfchofac (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), SF_NYCOEFF(sf),
YMATRIX(SF_YMATRIX(sf)), ier)
if (SF_XTERMS(sf) == YES) {
# solve for the nxcoeff right sides
cptr = SF_COEFF(sf)
do i = 1, SF_NXCOEFF(sf) {
call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
cptr = cptr + SF_NYCOEFF(sf)
}
} else {
cptr = SF_COEFF(sf)
call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
do i = 2, SF_NXCOEFF(sf) {
cptr = cptr + SF_NYCOEFF(sf)
COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf)
}
}
end
|
