# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include <math/surfit.h>
include "surfitdef.h"

# ISLFIT -- Procedure to fit a single line of a surface. The inner products
# of the x basis functions are calculated and accumulated into the
# SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The main diagonal is stored
# in the first row of XMATRIX followed by the non-zero lower diagonals. This
# method of storage is very efficient for the large symmetric banded matrices
# generated by spline fits. The Cholesky factorization of XMATRIX is calculated
# and stored in XMATRIX destroying the original data. The inner products
# of the x basis functions and the data ordinates are stored in the lineno-th
# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The x coefficients
# for each line are calculated by forward and back substitution and
# stored in the lineno-th line of XCOEFF destroying the original data.

procedure islfit (sf, cols, lineno, z, w, ncols, wtflag, ier)

pointer	sf		# pointer to the surface descriptor
int	cols[ncols]	# array of columns
int	lineno    	# lineno
real	z[ncols]	# the surface values
real	w[ncols]	# array of weights
int	ncols		# the number of columns
int	wtflag		# type of weighting
int	ier		# error codes

int	i, ii, j, k
pointer	xbzptr, xbptr
pointer	xmzptr, xmindex
pointer	xczptr, xcindex
pointer	xlzptr
pointer	left, rows, bw
pointer	sp
real	adotr()

begin
	# index the pointers
	xbzptr = SF_XBASIS(sf) - 1
	xmzptr = SF_XMATRIX(sf)
	xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1

	# zero the accumulators
	call aclrr (XMATRIX(SF_XMATRIX(sf)), SF_NXCOEFF(sf) * SF_XORDER(sf))
	call aclrr (XCOEFF(xczptr + 1), SF_NXCOEFF(sf))

	# free space used by previous islrefit calls
	if (SF_WZ(sf) != NULL)
	    call mfree (SF_WZ(sf), MEM_TYPE)
	if (SF_TLEFT(sf) != NULL)
	    call mfree (SF_TLEFT(sf), TY_INT)

	# reset number of points
	SF_NXPTS(sf) = ncols

	# calculate the weights, default is uniform weighting
	switch (wtflag) {
	case SF_UNIFORM:
	    call amovkr (1.0, w, ncols)
	case SF_USER:
	    # do not alter weights
	default:
	    call amovkr (1.0, w, ncols)
	}

	# allocate temporary space
	call smark (sp)
	call salloc (bw, ncols, TY_REAL)

	switch (SF_TYPE(sf)) {
	case SF_LEGENDRE, SF_CHEBYSHEV:

	    do i = 1, SF_XORDER(sf) {
		do j = 1, ncols
		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])

		xcindex = xczptr + i
		XCOEFF(xcindex) = XCOEFF(xcindex) + adotr (Memr[bw], z, ncols)

		xbptr = xbzptr
		ii = 0
		do k = i, SF_XORDER(sf) {
		    xmindex = xmzptr + ii
		    do j = 1, ncols
		        XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
			    XBASIS(xbptr+cols[j])
		    ii = ii + 1
		    xbptr = xbptr + SF_NCOLS(sf)
		}

		xbzptr = xbzptr + SF_NCOLS(sf)
		xmzptr = xmzptr + SF_XORDER(sf)
	    }

	case SF_SPLINE3, SF_SPLINE1:
	    xlzptr = SF_XLEFT(sf) - 1

	    call salloc (left, ncols, TY_INT)
	    call salloc (rows, ncols, TY_INT)

	    do j = 1, ncols
		Memi[left+j-1] = XLEFT(xlzptr+cols[j])
	    call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols)
	    call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols)
	    call aaddki (Memi[left], xczptr, Memi[left], ncols)

	    do i = 1, SF_XORDER(sf) {
		do j = 1, ncols {
		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
		    xcindex = Memi[left+j-1] + i
		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
		}

		xbptr = xbzptr
		ii = 0
		do k = i, SF_XORDER(sf) {
		    do j = 1, ncols {
			xmindex = Memi[rows+j-1] + ii
			XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
			    XBASIS(xbptr+cols[j])
		    }
		    ii = ii + 1
		    xbptr = xbptr + SF_NCOLS(sf)
		}

		xbzptr = xbzptr + SF_NCOLS(sf)
		call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols)
	    }
	}

	# release space
	call sfree (sp)

	# return if not enough data points
	ier = OK
	if ((SF_NXPTS(sf) - SF_NXCOEFF(sf)) < 0) {
	    ier = NO_DEG_FREEDOM
	    return
	}

	# calculate the Cholesky factorization of XMATRIX
	call sfchofac (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
	    XMATRIX(SF_XMATRIX(sf)), ier)

	# calculate the x coefficients for lineno-th image line and place in the
	# lineno-th row of XCOEFF
	call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
	    XCOEFF(xczptr + 1), XCOEFF(xczptr + 1))
end