1 files changed, 91 insertions, 0 deletions
diff --git a/math/gsurfit/gsevalr.x b/math/gsurfit/gsevalr.x
new file mode 100644
index 00000000..738e9915
--- /dev/null
+++ b/math/gsurfit/gsevalr.x
@@ -0,0 +1,91 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/gsurfit.h>
+include "gsurfitdef.h"
+
+# GSEVAL -- Procedure to evaluate the fitted surface at a single point.
+# The GS_NCOEFF(sf) coefficients are stored in the vector COEFF.
+
+real procedure gseval (sf, x, y)
+
+pointer	sf		# pointer to surface descriptor structure
+real	x		# x value
+real	y		# y value
+
+real	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_REAL)
+	    call salloc (yb, GS_NYCOEFF(sf), TY_REAL)
+	    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 rgs_b1cheb (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf),
+		XBS(xb))
+	    call rgs_b1cheb (y, GS_NYCOEFF(sf), GS_YMAXMIN(sf), GS_YRANGE(sf),
+		YBS(yb))
+	case GS_LEGENDRE:
+	    call rgs_b1leg (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf),
+		XBS(xb))
+	    call rgs_b1leg (y, GS_NYCOEFF(sf), GS_YMAXMIN(sf), GS_YRANGE(sf),
+		YBS(yb))
+	case GS_POLYNOMIAL:
+	    call rgs_b1pol (x, GS_NXCOEFF(sf), GS_XMAXMIN(sf), GS_XRANGE(sf),
+		XBS(xb))
+	    call rgs_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