1 files changed, 117 insertions, 0 deletions
diff --git a/math/gsurfit/gsfit1.gx b/math/gsurfit/gsfit1.gx
new file mode 100644
index 00000000..4e7341e4
--- /dev/null
+++ b/math/gsurfit/gsfit1.gx
@@ -0,0 +1,117 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/gsurfit.h>
+$if (datatype == r)
+include "gsurfitdef.h"
+$else
+include "dgsurfitdef.h"
+$endif
+
+# GSFIT1 -- Procedure to solve the normal equations for a surface.
+#
+# This version modifies the fitting matrix to remove the first
+# term from the fitting.  For the polynomial functions this means
+# constraining the constant term to be zero.  Note that the first
+# coefficent is still returned but with a value of zero.
+
+$if (datatype == r)
+procedure gsfit1 (sf, x, y, z, w, npts, wtflag, ier)
+$else
+procedure dgsfit1 (sf, x, y, z, w, npts, wtflag, ier)
+$endif
+
+pointer	sf		# surface descriptor
+PIXEL	x[npts]		# array of x values
+PIXEL	y[npts]		# array of y values
+PIXEL	z[npts]		# data array
+PIXEL	w[npts]		# array of weights
+int	npts		# number of data points
+int	wtflag		# type of weighting
+int	ier		# ier = OK, everything OK
+			# ier = SINGULAR, matrix is singular, 1 or more
+			# coefficients are 0.
+			# ier = NO_DEG_FREEDOM, too few points to solve matrix
+
+begin
+	$if (datatype == r)
+	    call gszero (sf)
+	    call gsacpts (sf, x, y, z, w, npts, wtflag)
+	    call gssolve1 (sf, ier)
+	$else
+	    call dgszero (sf)
+	    call dgsacpts (sf, x, y, z, w, npts, wtflag)
+	    call dgssolve1 (sf, ier)
+	$endif
+end
+
+
+# GSSOLVE1 -- Solve the matrix normal equations of the form ca = b for
+# a, where c is a symmetric, positive semi-definite, banded matrix with
+# GS_NXCOEFF(sf) * GS_NYCOEFF(sf) rows and a and b are GS_NXCOEFF(sf) *
+# GS_NYCOEFF(sf)-vectors.  Initially c is stored in the matrix MATRIX and b
+# is stored in VECTOR.  The Cholesky factorization of MATRIX is calculated
+# and stored in CHOFAC.  Finally the coefficients are calculated by forward
+# and back substitution and stored in COEFF.
+#
+# This version modifies the fitting matrix to remove the first
+# term from the fitting.  For the polynomial functions this means
+# constraining the constant term to be zero.  Note that the first
+# coefficent is still returned but with a value of zero.
+
+$if (datatype == r)
+procedure gssolve1 (sf, ier)
+$else
+procedure dgssolve1 (sf, ier)
+$endif
+
+pointer	sf 		# curve descriptor
+int	ier		# ier = OK, everything OK
+			# ier = SINGULAR, matrix is singular, 1 or more
+			# coefficients are 0.
+			# ier = NO_DEG_FREEDOM, too few points to solve matrix
+
+int	i, ncoeff, offset
+pointer	sp, vector, matrix
+
+begin
+
+	# test for number of degrees of freedom
+	offset = 1
+	ncoeff = GS_NCOEFF(sf) - offset
+	ier = OK
+	i = GS_NPTS(sf) - ncoeff
+	if (i < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# allocate working space for the reduced vector and matrix
+	call smark (sp)
+	call salloc (vector, ncoeff, TY_PIXEL)
+	call salloc (matrix, ncoeff*ncoeff, TY_PIXEL)
+
+	# eliminate the first term from the vector and matrix
+	call amov$t (VECTOR(GS_VECTOR(sf)+offset), Mem$t[vector], ncoeff)
+	do i = 0, ncoeff-1
+	    call amov$t (MATRIX(GS_MATRIX(sf)+(i+offset)*GS_NCOEFF(sf)),
+		Mem$t[matrix+i*ncoeff], ncoeff)
+
+	# solve for the coefficients.
+	switch (GS_TYPE(sf)) {
+	case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
+
+	    # calculate the Cholesky factorization of the data matrix
+	    call $tgschofac (Memd[matrix], ncoeff, ncoeff,
+	        CHOFAC(GS_CHOFAC(sf)), ier)
+
+	    # solve for the coefficients by forward and back substitution
+	    COEFF(GS_COEFF(sf)) = 0.
+	    call $tgschoslv (CHOFAC(GS_CHOFAC(sf)), ncoeff, ncoeff,
+	        Memd[vector], COEFF(GS_COEFF(sf)+offset))
+
+	default:
+	    call error (0, "GSSOLVE1: Illegal surface type.")
+	}
+
+	call sfree (sp)
+end