# 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