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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/surfit.h>
include "surfitdef.h"
# ISLREFIT -- Procedure to refit the data assuming that the cols and w
# arrays do not change. SIFLREFIT assumes that the Cholesky factorization
# of XMATRIX is stored in XMATRIX. The inner products of the x basis
# functions and the data ordinates are accumulated into the lineno-th
# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The coefficients
# for line number lineno are calculated by forward and back
# substitution and placed in the lineno-th row of XCOEFF replacing the
# original data.
procedure islrefit (sf, cols, lineno, z, w)
pointer sf # pointer to surface descriptor
int cols[ARB] # columns to be fit
int lineno # line number
real z[ARB] # surface values
real w[ARB] # weight values
int i, j
pointer xbzptr, xczptr, xcindex, xlzptr
begin
# set pointers
xbzptr = SF_XBASIS(sf) - 1
xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
xlzptr = SF_XLEFT(sf) - 1
# reset lineno-th row of the x coefficient matrix
call aclrr (XCOEFF(xczptr+1), SF_NXCOEFF(sf))
if (SF_WZ(sf) == NULL)
call malloc (SF_WZ(sf), SF_NXPTS(sf), MEM_TYPE)
# calculate new right sides
call amulr (w, z, Memr[SF_WZ(sf)], SF_NXPTS(sf))
switch (SF_TYPE(sf)) {
case SF_LEGENDRE, SF_CHEBYSHEV:
do i = 1, SF_XORDER(sf) {
xcindex = xczptr + i
do j = 1, SF_NXPTS(sf)
XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] *
XBASIS(xbzptr+cols[j])
xbzptr = xbzptr + SF_NCOLS(sf)
}
case SF_SPLINE3, SF_SPLINE1:
if (SF_TLEFT(sf) == NULL)
call malloc (SF_TLEFT(sf), SF_NXPTS(sf), TY_INT)
do i = 1, SF_NXPTS(sf)
Memi[SF_TLEFT(sf)+i-1] = XLEFT(xlzptr+cols[i]) + xczptr
do i = 1, SF_XORDER(sf) {
do j = 1, SF_NXPTS(sf) {
xcindex = Memi[SF_TLEFT(sf)+j-1] + i
XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] *
XBASIS(xbzptr+cols[j])
}
xbzptr = xbzptr + SF_NCOLS(sf)
}
}
# solve for the new x coefficients for line number lineno
call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
XCOEFF(xczptr+1), XCOEFF(xczptr+1))
end
|
