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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/surfit.h>
include "surfitdef.h"
# ISLFIT -- Procedure to fit a single line of a surface. The inner products
# of the x basis functions are calculated and accumulated into the
# SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The main diagonal is stored
# in the first row of XMATRIX followed by the non-zero lower diagonals. This
# method of storage is very efficient for the large symmetric banded matrices
# generated by spline fits. The Cholesky factorization of XMATRIX is calculated
# and stored in XMATRIX destroying the original data. The inner products
# of the x basis functions and the data ordinates are stored in the lineno-th
# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The x coefficients
# for each line are calculated by forward and back substitution and
# stored in the lineno-th line of XCOEFF destroying the original data.
procedure islfit (sf, cols, lineno, z, w, ncols, wtflag, ier)
pointer sf # pointer to the surface descriptor
int cols[ncols] # array of columns
int lineno # lineno
real z[ncols] # the surface values
real w[ncols] # array of weights
int ncols # the number of columns
int wtflag # type of weighting
int ier # error codes
int i, ii, j, k
pointer xbzptr, xbptr
pointer xmzptr, xmindex
pointer xczptr, xcindex
pointer xlzptr
pointer left, rows, bw
pointer sp
real adotr()
begin
# index the pointers
xbzptr = SF_XBASIS(sf) - 1
xmzptr = SF_XMATRIX(sf)
xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
# zero the accumulators
call aclrr (XMATRIX(SF_XMATRIX(sf)), SF_NXCOEFF(sf) * SF_XORDER(sf))
call aclrr (XCOEFF(xczptr + 1), SF_NXCOEFF(sf))
# free space used by previous islrefit calls
if (SF_WZ(sf) != NULL)
call mfree (SF_WZ(sf), MEM_TYPE)
if (SF_TLEFT(sf) != NULL)
call mfree (SF_TLEFT(sf), TY_INT)
# reset number of points
SF_NXPTS(sf) = ncols
# calculate the weights, default is uniform weighting
switch (wtflag) {
case SF_UNIFORM:
call amovkr (1.0, w, ncols)
case SF_USER:
# do not alter weights
default:
call amovkr (1.0, w, ncols)
}
# allocate temporary space
call smark (sp)
call salloc (bw, ncols, TY_REAL)
switch (SF_TYPE(sf)) {
case SF_LEGENDRE, SF_CHEBYSHEV:
do i = 1, SF_XORDER(sf) {
do j = 1, ncols
Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
xcindex = xczptr + i
XCOEFF(xcindex) = XCOEFF(xcindex) + adotr (Memr[bw], z, ncols)
xbptr = xbzptr
ii = 0
do k = i, SF_XORDER(sf) {
xmindex = xmzptr + ii
do j = 1, ncols
XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
XBASIS(xbptr+cols[j])
ii = ii + 1
xbptr = xbptr + SF_NCOLS(sf)
}
xbzptr = xbzptr + SF_NCOLS(sf)
xmzptr = xmzptr + SF_XORDER(sf)
}
case SF_SPLINE3, SF_SPLINE1:
xlzptr = SF_XLEFT(sf) - 1
call salloc (left, ncols, TY_INT)
call salloc (rows, ncols, TY_INT)
do j = 1, ncols
Memi[left+j-1] = XLEFT(xlzptr+cols[j])
call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols)
call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols)
call aaddki (Memi[left], xczptr, Memi[left], ncols)
do i = 1, SF_XORDER(sf) {
do j = 1, ncols {
Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
xcindex = Memi[left+j-1] + i
XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
}
xbptr = xbzptr
ii = 0
do k = i, SF_XORDER(sf) {
do j = 1, ncols {
xmindex = Memi[rows+j-1] + ii
XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
XBASIS(xbptr+cols[j])
}
ii = ii + 1
xbptr = xbptr + SF_NCOLS(sf)
}
xbzptr = xbzptr + SF_NCOLS(sf)
call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols)
}
}
# release space
call sfree (sp)
# return if not enough data points
ier = OK
if ((SF_NXPTS(sf) - SF_NXCOEFF(sf)) < 0) {
ier = NO_DEG_FREEDOM
return
}
# calculate the Cholesky factorization of XMATRIX
call sfchofac (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
XMATRIX(SF_XMATRIX(sf)), ier)
# calculate the x coefficients for lineno-th image line and place in the
# lineno-th row of XCOEFF
call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
XCOEFF(xczptr + 1), XCOEFF(xczptr + 1))
end
|
