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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/surfit.h>
include "surfitdef.h"
# ISRESOLVE -- Procedure to solve the x coefficient matrix for the surface
# coefficients assuming that the lines array is unchanged since the last
# call to SIFSOLVE. The Cholesky factorization of YMATRIX is assumed to be
# stored in YMATRIX. The inner product of the y basis functions and i-th
# column of XCOEFF containing the i-th x coefficients for each line are
# calculated and stored in the i-th row of the SF_NYCOEFF(sf) by
# SF_NXCOEFF(sf) array COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF
# is solved to determine the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) surface
# coefficients. After a call to SIFSOLVE the coefficient for the i-th
# y basis function and the j-th x coefficient is found in the j-th row
# and i-th column of COEFF.
procedure isresolve (sf, lines, ier)
pointer sf # pointer to the surface descriptor structure
int lines[ARB] # line numbers included in the fit
int ier # error code
int i, j, k, nxcoeff
pointer ybzptr
pointer ylzptr
pointer xczptr, xcptr, xcindex
pointer czptr, cptr
pointer left, tleft
pointer sp
begin
# define pointers
ybzptr = SF_YBASIS(sf) - 1
xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1
czptr = SF_COEFF(sf) - 1
# zero out coefficient matrix
call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf))
switch (SF_TYPE(sf)) {
case SF_LEGENDRE, SF_CHEBYSHEV:
# loop over the y basis functions
nxcoeff = SF_NXCOEFF(sf)
do i = 1, SF_YORDER(sf) {
cptr = czptr + i
do k = 1, nxcoeff {
xcptr = xczptr + k
do j = 1, SF_NYPTS(sf) {
xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
COEFF(cptr) = COEFF(cptr) +
YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex)
}
cptr = cptr + SF_NYCOEFF(sf)
}
ybzptr = ybzptr + SF_NYPTS(sf)
# if SF_XTERMS(sf) = NO do not accumulate elements
# of COEFF where i != 1 and k != 1
if (SF_XTERMS(sf) == NO)
nxcoeff = 1
}
case SF_SPLINE3, SF_SPLINE1:
call smark (sp)
call salloc (left, SF_NYPTS(sf), TY_INT)
call salloc (tleft, SF_NYPTS(sf), TY_INT)
ylzptr = SF_YLEFT(sf) - 1
do j = 1, SF_NYPTS(sf)
Memi[left+j-1] = YLEFT(ylzptr+lines[j])
call aaddki (Memi[left], czptr, Memi[left], SF_NYPTS(sf))
nxcoeff = SF_NXCOEFF(sf)
do i = 1, SF_YORDER(sf) {
call aaddki (Memi[left], i, Memi[tleft], SF_NYPTS(sf))
do k = 1, nxcoeff {
xcptr = xczptr + k
do j = 1, SF_NYPTS(sf) {
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],
SF_NYPTS(sf))
}
ybzptr = ybzptr + SF_NYPTS(sf)
}
call sfree (sp)
}
# return if not enough points
ier = OK
if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) {
ier = NO_DEG_FREEDOM
return
}
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 {
# solve for the y coefficients
cptr = SF_NYCOEFF(sf)
call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
# solve for the x coefficients
do i = 2, SF_NXCOEFF(sf) {
cptr = czptr + SF_NYCOEFF(sf)
COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf)
}
}
end
|
