diff options
author | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
---|---|---|
committer | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
commit | fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch) | |
tree | bdda434976bc09c864f2e4fa6f16ba1952b1e555 /math/gsurfit/gseval.gx | |
download | iraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz |
Initial commit
Diffstat (limited to 'math/gsurfit/gseval.gx')
-rw-r--r-- | math/gsurfit/gseval.gx | 104 |
1 files changed, 104 insertions, 0 deletions
diff --git a/math/gsurfit/gseval.gx b/math/gsurfit/gseval.gx new file mode 100644 index 00000000..d57d21c7 --- /dev/null +++ b/math/gsurfit/gseval.gx @@ -0,0 +1,104 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/gsurfit.h> +$if (datatype == r) +include "gsurfitdef.h" +$else +include "dgsurfitdef.h" +$endif + +# GSEVAL -- Procedure to evaluate the fitted surface at a single point. +# The GS_NCOEFF(sf) coefficients are stored in the vector COEFF. + +$if (datatype == r) +real procedure gseval (sf, x, y) +$else +double procedure dgseval (sf, x, y) +$endif + +pointer sf # pointer to surface descriptor structure +PIXEL x # x value +PIXEL y # y value + +PIXEL 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: + $if (datatype == r) + call salloc (xb, GS_NXCOEFF(sf), TY_REAL) + call salloc (yb, GS_NYCOEFF(sf), TY_REAL) + $else + call salloc (xb, GS_NXCOEFF(sf), TY_DOUBLE) + call salloc (yb, GS_NYCOEFF(sf), TY_DOUBLE) + $endif + 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 $tgs_b1cheb (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf), + XBS(xb)) + call $tgs_b1cheb (y, GS_NYCOEFF(sf), GS_YMAXMIN(sf), GS_YRANGE(sf), + YBS(yb)) + case GS_LEGENDRE: + call $tgs_b1leg (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf), + XBS(xb)) + call $tgs_b1leg (y, GS_NYCOEFF(sf), GS_YMAXMIN(sf), GS_YRANGE(sf), + YBS(yb)) + case GS_POLYNOMIAL: + call $tgs_b1pol (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf), + XBS(xb)) + call $tgs_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 |