Initial commit

1 files changed, 147 insertions, 0 deletions

diff --git a/math/curfit/cv_beval.gx b/math/curfit/cv_beval.gx
new file mode 100644
index 00000000..3bb2b2e1
--- /dev/null
+++ b/math/curfit/cv_beval.gx
@@ -0,0 +1,147 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# CV_BCHEB -- Procedure to evaluate all the non-zero Chebyshev functions for
+# a set of points and given order.
+
+procedure $tcv_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)
+		call amovk$t (PIXEL(1.0), basis, npts)
+	    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)
+		call amulk$t (basis[bptr], PIXEL(2.0), basis[bptr], npts)
+		call asub$t (basis[bptr], basis[bptr-2*npts], basis[bptr], npts)
+	    }
+	    bptr = bptr + npts
+	}
+end
+
+
+# CV_BLEG -- Procedure to evaluate all the non zero Legendre function
+# for a given order and set of points.
+
+procedure $tcv_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)
+		call amovk$t (PIXEL(1.0), basis, npts)
+	    else if (k == 2)
+		call alta$t (x, basis[bptr], npts, k1, k2)
+	    else {
+		ri = k
+		ri1 = (PIXEL(2.0) * ri - PIXEL(3.0)) / (ri - PIXEL(1.0))
+		ri2 = - (ri - PIXEL(2.0)) / (ri - PIXEL(1.0))
+		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
+
+
+# CV_BSPLINE1 -- Evaluate all the non-zero spline1 functions for a set
+# of points.
+
+procedure $tcv_bspline1 (x, npts, npieces, k1, k2, basis, left)
+
+PIXEL	x[npts]		# set of data points
+int	npts		# number of points
+int	npieces		# number of polynomial pieces minus 1
+PIXEL	k1, k2		# normalizing constants
+PIXEL	basis[ARB]	# basis functions
+int	left[ARB]	# indices of the appropriate spline functions
+
+int	k
+
+begin
+	call alta$t (x, basis[1+npts], npts, k1, k2)
+	call acht$ti (basis[1+npts], left, npts)
+	call aminki (left, npieces, left, npts)
+
+	do k = 1, npts {
+	    basis[npts+k] = max (PIXEL(0.0), min (PIXEL(1.0),
+	        basis[npts+k] - left[k]))
+	    basis[k] = max (PIXEL(0.0), min (PIXEL(1.0), PIXEL(1.0) -
+	        basis[npts+k]))
+	}
+end
+
+
+# CV_BSPLINE3 --  Procedure to evaluate all the non-zero basis functions
+# for a cubic spline.
+
+procedure $tcv_bspline3 (x, npts, npieces, k1, k2, basis, left)
+
+PIXEL	x[npts]		# array of data points
+int	npts		# number of data points
+int	npieces		# number of polynomial pieces minus 1
+PIXEL	k1, k2		# normalizing constants
+PIXEL	basis[ARB]	# array of basis functions
+int	left[ARB]	# array of indices for first non-zero spline
+
+int	i
+pointer	sp, sx, tx
+PIXEL	dsx, dtx
+
+begin
+	# allocate space
+	call smark (sp)
+	call salloc (sx, npts, TY_PIXEL)
+	call salloc (tx, npts, TY_PIXEL)
+
+	# calculate the index of the first non-zero coeff
+	call alta$t (x, Mem$t[sx], npts, k1, k2)
+	call acht$ti (Mem$t[sx], left, npts)
+	call aminki (left, npieces, left, npts)
+
+	do i = 1, npts {
+	    Mem$t[sx+i-1] = max (PIXEL(0.0), min (PIXEL(1.0),
+	        Mem$t[sx+i-1] - left[i]))
+	    Mem$t[tx+i-1] = max (PIXEL(0.0), min (PIXEL(1.0), PIXEL(1.0) -
+	        Mem$t[sx+i-1]))
+	}
+
+	# calculate the basis function
+	#call apowk$t (Mem$t[tx], 3, basis, npts)
+	do i = 1, npts {
+	    dsx = Mem$t[sx+i-1]
+	    dtx = Mem$t[tx+i-1]
+	    basis[i] = dtx * dtx * dtx
+	    basis[npts+i] = PIXEL(1.0) + dtx * (PIXEL(3.0) + dtx *
+	        (PIXEL(3.0) - PIXEL(3.0) * dtx))
+	    basis[2*npts+i] = PIXEL(1.0) + dsx * (PIXEL(3.0) + dsx *
+	        (PIXEL(3.0) - PIXEL(3.0) * dsx))
+	    basis[3*npts+i] = dsx * dsx * dsx
+	}
+	#call apowk$t (Mem$t[sx], 3, basis[1+3*npts], npts)
+
+	# release space
+	call sfree (sp)
+end