1 files changed, 120 insertions, 0 deletions
diff --git a/math/gsurfit/gs_beval.gx b/math/gsurfit/gs_beval.gx
new file mode 100644
index 00000000..da45f122
--- /dev/null
+++ b/math/gsurfit/gs_beval.gx
@@ -0,0 +1,120 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# GS_BPOL -- Procedure to evaluate all the non-zero polynomial functions for
+# a set of points and given order.
+
+procedure $tgs_bpol (x, npts, order, k1, k2, basis)
+
+PIXEL	x[npts]		# array of data points
+int	npts		# number of points
+int	order		# order of polynomial, order = 1, constant
+PIXEL	k1, k2		# normalizing constants
+PIXEL	basis[ARB]	# basis functions
+
+int	bptr, k
+
+begin
+	bptr = 1
+	do k = 1, order {
+
+	    if (k == 1)
+		$if (datatype == r)
+	        call amovkr (1.0, basis, npts)
+		$else
+	        call amovkd (1.0d0, basis, npts)
+		$endif
+	    else if (k == 2)
+		call amov$t (x, basis[bptr], npts)
+	    else 
+		call amul$t (basis[bptr-npts], x, basis[bptr], npts)
+		
+	    bptr = bptr + npts
+	}
+end
+
+# GS_BCHEB -- Procedure to evaluate all the non-zero Chebyshev functions for
+# a set of points and given order.
+
+procedure $tgs_bcheb (x, npts, order, k1, k2, basis)
+
+PIXEL	x[npts]		# array of data points
+int	npts		# number of points
+int	order		# order of polynomial, order = 1, constant
+PIXEL	k1, k2		# normalizing constants
+PIXEL	basis[ARB]	# basis functions
+
+int	k, bptr
+
+begin
+	bptr = 1
+	do k = 1, order {
+
+	    if (k == 1)
+		$if (datatype == r)
+	        call amovkr (1.0, basis, npts)
+		$else
+		call amovkd (1.0d0, basis, npts)
+		$endif
+	    else if (k == 2)
+		call alta$t (x, basis[bptr], npts, k1, k2)
+	    else {
+		call amul$t (basis[1+npts], basis[bptr-npts], basis[bptr],
+		    npts)
+		$if (datatype == r)
+		call amulkr (basis[bptr], 2.0, basis[bptr], npts)
+		$else
+		call amulkd (basis[bptr], 2.0d0, basis[bptr], npts)
+		$endif
+		call asub$t (basis[bptr], basis[bptr-2*npts], basis[bptr], npts)
+	    }
+		
+	    bptr = bptr + npts
+	}
+end
+
+
+# GS_BLEG -- Procedure to evaluate all the non zero Legendre function
+# for a given order and set of points.
+
+procedure $tgs_bleg (x, npts, order, k1, k2, basis)
+
+PIXEL	x[npts]		# number of data points
+int	npts		# number of points
+int	order		# order of polynomial, 1 is a constant
+PIXEL	k1, k2		# normalizing constants
+PIXEL	basis[ARB]	# array of basis functions
+
+int	k, bptr
+PIXEL	ri, ri1, ri2
+
+begin
+	bptr = 1
+	do k = 1, order {
+
+	    if (k == 1)
+		$if (datatype == r)
+		call amovkr (1.0, basis, npts)
+		$else
+		call amovkd (1.0d0, basis, npts)
+		$endif
+	    else if (k == 2)
+		call alta$t (x, basis[bptr], npts, k1, k2)
+	    else {
+		$if (datatype == r)
+		ri = k
+		ri1 = (2. * ri - 3.) / (ri - 1.)
+		ri2 = - (ri - 2.) / (ri - 1.)
+		$else
+		ri = k
+		ri1 = (2.0d0 * ri - 3.0d0) / (ri - 1.0d0)
+		ri2 = - (ri - 2.0d0) / (ri - 1.0d0)
+		$endif
+		call amul$t (basis[1+npts], basis[bptr-npts], basis[bptr],
+		    npts)
+		call awsu$t (basis[bptr], basis[bptr-2*npts],
+			basis[bptr], npts, ri1, ri2)
+	    }
+			
+	    bptr = bptr + npts
+	}
+end