Repatch (from linux) of OSX IRAF

1 files changed, 174 insertions, 0 deletions

diff --git a/math/gsurfit/gsrefit.gx b/math/gsurfit/gsrefit.gx
new file mode 100644
index 00000000..00327abb
--- /dev/null
+++ b/math/gsurfit/gsrefit.gx
@@ -0,0 +1,174 @@
+# 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
+
+# GSREFIT -- Procedure to refit the surface assuming that the x, y and w
+# values and the matrices MATRIX and CHOFAC have remained unchanged. It
+# is necessary only to accumulate a new VECTOR. The new coefficients
+# are calculated by forward and back subsitution and stored in COEFF.
+
+$if (datatype == r)
+procedure gsrefit (sf, x, y, z, w, ier)
+$else
+procedure dgsrefit (sf, x, y, z, w, ier)
+$endif
+
+pointer	sf		# surface descriptor
+PIXEL	x[ARB]		# array of x values
+PIXEL	y[ARB]		# array of y values
+PIXEL	z[ARB]		# data array
+PIXEL	w[ARB]		# array of weights
+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	k, l
+int	xorder, nfree, maxorder
+pointer	sp, vzptr, vindex, bxptr, byptr, bwz
+
+PIXEL	adot$t()
+
+errchk	smark, salloc, sfree
+
+begin
+	# clear accumulator
+	call aclr$t (VECTOR(GS_VECTOR(sf)), GS_NCOEFF(sf))
+
+	# if first call to gsefit calculate basis functions
+	if (GS_XBASIS(sf) == NULL || GS_YBASIS(sf) == NULL) {
+
+	    $if (datatype == r)
+	    call malloc (GS_WZ(sf), GS_NPTS(sf), TY_REAL)
+	    $else
+	    call malloc (GS_WZ(sf), GS_NPTS(sf), TY_DOUBLE)
+	    $endif
+
+	    switch (GS_TYPE(sf)) {
+	    case GS_LEGENDRE:
+		$if (datatype == r)
+	        call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
+		    TY_REAL)
+	        call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
+		    TY_REAL)
+		$else
+	        call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
+		    TY_DOUBLE)
+	        call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
+		    TY_DOUBLE)
+		$endif
+	        call $tgs_bleg (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
+	    		  GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
+	        call $tgs_bleg (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
+	    		  GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
+	    case GS_CHEBYSHEV:
+		$if (datatype == r)
+	        call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
+		    TY_REAL)
+	        call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
+		    TY_REAL)
+		$else
+	        call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
+		    TY_DOUBLE)
+	        call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
+		    TY_DOUBLE)
+		$endif
+	        call $tgs_bcheb (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
+	    		  GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
+	        call $tgs_bcheb (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
+	    		  GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
+	    case GS_POLYNOMIAL:
+		$if (datatype == r)
+	        call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
+		    TY_REAL)
+	        call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
+		    TY_REAL)
+		$else
+	        call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
+		    TY_DOUBLE)
+	        call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
+		    TY_DOUBLE)
+		$endif
+	        call $tgs_bpol (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
+	    		  GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
+	        call $tgs_bpol (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
+	    		  GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
+	    default:
+		call error (0, "GSREFIT: Unknown curve type.")
+	    }
+
+	}
+
+	call smark (sp)
+	$if (datatype == r)
+	call salloc (bwz, GS_NPTS(sf), TY_REAL) 
+	$else
+	call salloc (bwz, GS_NPTS(sf), TY_DOUBLE)
+	$endif
+
+	# index the pointers
+	vzptr = GS_VECTOR(sf) - 1
+	byptr = GS_YBASIS(sf)
+
+	switch (GS_TYPE(sf)) {
+	case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
+
+	    call amul$t (w, z, Mem$t[GS_WZ(sf)], GS_NPTS(sf))
+	    xorder = GS_XORDER(sf)
+	    maxorder = max (GS_XORDER(sf) + 1, GS_YORDER(sf) + 1)
+
+	    do l = 1, GS_YORDER(sf) {
+		call amul$t (Mem$t[GS_WZ(sf)], YBASIS(byptr), Mem$t[bwz],
+		    GS_NPTS(sf))
+		bxptr = GS_XBASIS(sf)
+	        do k = 1, xorder {
+		    vindex = vzptr + k
+	            VECTOR(vindex) = VECTOR(vindex) + adot$t (Mem$t[bwz],
+		        XBASIS(bxptr), GS_NPTS(sf))
+		    bxptr = bxptr + GS_NPTS(sf)
+	        }
+
+		vzptr = vzptr + xorder
+		switch (GS_XTERMS(sf)) {
+		case GS_XNONE:
+		    xorder = 1
+		case GS_XHALF:
+		    if ((l + GS_XORDER(sf) + 1) > maxorder)
+		        xorder = xorder - 1
+		default:
+		    ;
+		}
+		byptr = byptr + GS_NPTS(sf)
+	    }
+
+	default:
+	    call error (0, "GSACCUM: Unknown curve type.")
+	}
+
+	# test for number of degrees of freedom
+	ier = OK
+	nfree = GS_NPTS(sf) - GS_NCOEFF(sf)
+	if (nfree < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# calculate the values of the coefficients
+	switch (GS_TYPE(sf)) {
+	case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
+
+	    # solve for the coefficients by forward and back substitution
+	    call $tgschoslv (CHOFAC(GS_CHOFAC(sf)), GS_NCOEFF(sf),
+	        GS_NCOEFF(sf), VECTOR(GS_VECTOR(sf)), COEFF(GS_COEFF(sf)))
+	default:
+	    call error (0, "GSSOLVE: Illegal surface type.")
+	}
+
+	# release the space
+	call sfree (sp)
+end