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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/gsurfit.h>
include "dgsurfitdef.h"
# GSEVAL -- Procedure to evaluate the fitted surface at a single point.
# The GS_NCOEFF(sf) coefficients are stored in the vector COEFF.
double procedure dgseval (sf, x, y)
pointer sf # pointer to surface descriptor structure
double x # x value
double y # y value
double sum, accum
int i, ii, k, maxorder, xorder
pointer sp, xb, xzb, yb, yzb, czptr
errchk smark, salloc, sfree
begin
call smark (sp)
# allocate space for the basis functions
switch (GS_TYPE(sf)) {
case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
call salloc (xb, GS_NXCOEFF(sf), TY_DOUBLE)
call salloc (yb, GS_NYCOEFF(sf), TY_DOUBLE)
xzb = xb - 1
yzb = yb - 1
czptr = GS_COEFF(sf) - 1
default:
call error (0, "GSEVAL: Unknown curve type.")
}
# calculate the basis functions
switch (GS_TYPE(sf)) {
case GS_CHEBYSHEV:
call dgs_b1cheb (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf),
XBS(xb))
call dgs_b1cheb (y, GS_NYCOEFF(sf), GS_YMAXMIN(sf), GS_YRANGE(sf),
YBS(yb))
case GS_LEGENDRE:
call dgs_b1leg (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf),
XBS(xb))
call dgs_b1leg (y, GS_NYCOEFF(sf), GS_YMAXMIN(sf), GS_YRANGE(sf),
YBS(yb))
case GS_POLYNOMIAL:
call dgs_b1pol (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf),
XBS(xb))
call dgs_b1pol (y, GS_NYCOEFF(sf), GS_YMAXMIN(sf), GS_YRANGE(sf),
YBS(yb))
default:
call error (0, "GSEVAL: Unknown surface type.")
}
# initialize accumulator
# basis functions
sum = 0.
# loop over y basis functions
maxorder = max (GS_XORDER(sf) + 1, GS_YORDER(sf) + 1)
xorder = GS_XORDER(sf)
ii = 1
do i = 1, GS_YORDER(sf) {
# loop over the x basis functions
accum = 0.
do k = 1, xorder {
accum = accum + COEFF(czptr+ii) * XBS(xzb+k)
ii = ii + 1
}
accum = accum * YBS(yzb+i)
sum = sum + accum
# elements of COEFF where neither k = 1 or i = 1
# are not calculated if GS_XTERMS(sf) = NO
switch (GS_XTERMS(sf)) {
case GS_XNONE:
xorder = 1
case GS_XHALF:
if ((i + GS_XORDER(sf) + 1) > maxorder)
xorder = xorder - 1
default:
;
}
}
call sfree (sp)
return (sum)
end
|
