Repatch (from linux) of OSX IRAF

1 files changed, 348 insertions, 0 deletions

diff --git a/math/gsurfit/zzdebug.x b/math/gsurfit/zzdebug.x
new file mode 100644
index 00000000..522da747
--- /dev/null
+++ b/math/gsurfit/zzdebug.x
@@ -0,0 +1,348 @@
+task test = t_test
+
+include <math/gsurfit.h>
+
+procedure t_test()
+
+int	i, j, k
+int	xorder, yorder, xterms, stype
+int	ncoeff, maxorder, xincr, npts, stype1, ier
+pointer	gs, ags, sgs
+double	dx, dy, const, accum, sum, rms1, rms2
+double	x[121], y[121], z[121], w[121], zfit[121], coeff[121], save[121]
+int	clgeti()
+double	dgseval(), dgsgcoeff()
+
+begin
+	# Generate x and y grid.
+	dy = -1.0d0
+	npts = 0
+	do i = 1, 9 {
+	    dx = -1.0d0
+	    do j = 1, 9 {
+		x[npts+1] = dx
+		y[npts+1] = dy
+		npts = npts + 1
+		dx = dx + 0.25d0
+	    }
+	    dy = dy + 0.25d0
+	}
+
+	stype = clgeti ("stype")
+	xorder = clgeti ("xorder")
+	yorder = clgeti ("yorder")
+	xterms = clgeti ("xterms")
+	call printf ("\n\nSURFACE: %d XORDER: %d YORDER: %d XTERMS: %d\n")
+	    call pargi (stype)
+	    call pargi (xorder)
+	    call pargi (yorder)
+	    call pargi (xterms)
+
+	# Generate data
+	if (stype > 3) {
+	    switch (xterms) {
+	    case GS_XNONE:
+
+		do i = 1, npts {
+		    sum = 0.0d0
+		    do j = 2, yorder
+			sum = sum + j * y[i] ** (j - 1)
+		    do j = 2, xorder
+			sum = sum + j * x[i] ** (j - 1)
+		    z[i] = sum
+		}
+
+	    case GS_XHALF:
+
+		maxorder = max (xorder + 1, yorder + 1)
+		do i = 1, npts {
+		    sum = 0.0d0
+		    xincr = xorder
+		    do j = 1, yorder {
+			const = j * y[i] ** (j - 1)
+			accum= 0.0d0
+			do k = 1, xincr {
+			    if (j > 1 || k > 1)
+				accum = accum + k * x[i] ** (k - 1)
+			}
+			sum = sum + const * accum
+			if ((j + xorder + 1) > maxorder)
+			    xincr = xincr - 1
+		    }
+		    z[i] = sum
+		}
+
+	    case GS_XFULL:
+
+		do i = 1, npts {
+		    sum = 0.0d0
+		    do j = 1, yorder {
+			const = j * y[i] ** (j - 1)
+			accum = 0.0d0
+			do k = 1, xorder {
+			    if (j > 1 || k > 1)
+				accum = accum + k * x[i] ** (k - 1)
+			}
+			sum = sum + const * accum
+		    }
+		    z[i] = sum
+		}
+	    }
+
+	    stype1 = stype - 3
+	} else {
+	    switch (xterms) {
+	    case GS_XNONE:
+
+		do i = 1, npts {
+		    sum = 0.0d0
+		    do j = 2, yorder
+			sum = sum + j * y[i] ** (j - 1)
+		    do j = 1, xorder
+			sum = sum + j * x[i] ** (j - 1)
+		    z[i] = sum
+		}
+
+	    case GS_XHALF:
+
+		maxorder = max (xorder + 1, yorder + 1)
+		do i = 1, npts {
+		    sum = 0.0d0
+		    xincr = xorder
+		    do j = 1, yorder {
+			const = j * y[i] ** (j - 1)
+			accum= 0.0d0
+			do k = 1, xincr {
+			    accum = accum + k * x[i] ** (k - 1)
+			}
+			sum = sum + const * accum
+			if ((j + xorder + 1) > maxorder)
+			    xincr = xincr - 1
+		    }
+		    z[i] = sum
+		}
+
+	    case GS_XFULL:
+
+		do i = 1, npts {
+		    sum = 0.0d0
+		    do j = 1, yorder {
+			const = j * y[i] ** (j - 1)
+			accum = 0.0d0
+			do k = 1, xorder {
+			    accum = accum + k * x[i] ** (k - 1)
+			}
+			sum = sum + const * accum
+		    }
+		    z[i] = sum
+		}
+	    }
+
+	    stype1 = stype
+	}
+
+	# Print out the data.
+	call printf ("\nXIN:\n")
+	do i = 1, npts {
+	    call printf ("%6.3f ")
+		call pargd (x[i])
+	    if (mod (i, 9) == 0)
+		call printf ("\n")
+	}
+	call printf ("\n")
+
+	call printf ("\nYIN:\n")
+	do i = 1, npts {
+	    call printf ("%6.3f ")
+		call pargd (y[i])
+	    if (mod (i, 9) == 0)
+		call printf ("\n")
+	}
+	call printf ("\n")
+
+	call printf ("\nZIN:\n")
+	do i = 1, npts {
+	    call printf ("%6.3f ")
+		call pargd (z[i])
+	    if (mod (i, 9) == 0)
+		call printf ("\n")
+	}
+	call printf ("\n")
+
+	# Fit surface.
+	call dgsinit (gs, stype1, xorder, yorder, xterms, -1.0d0, 1.0d0,
+	    -1.0d0, 1.0d0)
+	if (stype > 3) {
+	    call dgsset (gs, GSXREF, 0d0)
+	    call dgsset (gs, GSYREF, 0d0)
+	    call dgsset (gs, GSZREF, 0d0)
+	}
+	call dgsfit (gs, x, y, z, w, npts, WTS_UNIFORM, ier) 
+	call printf ("\nFIT ERROR CODE: %d\n")
+	    call pargi (ier)
+
+	# Evaluate the fit and its rms.
+	call dgsvector (gs, x, y, zfit, npts)
+	call printf ("\nZFIT:\n")
+	do i = 1, npts {
+	    call printf ("%6.3f ")
+		call pargd (zfit[i])
+	    if (mod (i, 9) == 0)
+		call printf ("\n")
+	}
+	call printf ("\n")
+	rms1 = 0.0d0
+	do i = 1, npts
+	    rms1 = rms1 + (z[i] - zfit[i]) ** 2 
+	rms1 = sqrt (rms1 / (npts - 1))
+	rms2 = 0.0d0
+	do i = 1, npts
+	    rms2 = rms2 + (z[i] - dgseval (gs, x[i], y[i])) ** 2
+	rms2 = sqrt (rms2 / (npts - 1))
+	#call printf ("\nRMS: vector = %0.14g point = %0.14g\n\n")
+	call printf ("\nRMS: vector = %0.4f point = %0.4f\n\n")
+	    call pargd (rms1)
+	    call pargd (rms2)
+
+	# Print the coefficients.
+	call dgscoeff (gs, coeff, ncoeff)
+	call printf ("GSFIT coeff:\n")
+	call printf ("first %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (gs, 1, 1))
+	    call pargd (dgsgcoeff (gs, xorder, 1))
+	do i = 1, ncoeff {
+	    call printf ("%d  %0.14g\n")
+		call pargi (i)
+		call pargd (coeff[i])
+	}
+	call printf ("last %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (gs, 1, yorder))
+	    call pargd (dgsgcoeff (gs, xorder, yorder))
+	call printf ("\n")
+
+	call dgsfree (gs)
+	return
+
+	# Evaluate the first derivatives.
+	call dgsder (gs, x, y, zfit, npts, 1, 0)
+	call printf ("\nZDER: 1 0\n")
+	do i = 1, npts {
+	    call printf ("%0.7g ")
+		call pargd (zfit[i])
+	    if (mod (i, 9) == 0)
+		call printf ("\n")
+	}
+	call printf ("\n")
+
+	call dgsder (gs, x, y, zfit, npts, 0, 1)
+	call printf ("\nZDER: 0 1\n")
+	do i = 1, npts {
+	    call printf ("%0.7g ")
+		call pargd (zfit[i])
+	    if (mod (i, 9) == 0)
+		call printf ("\n")
+	}
+	call printf ("\n")
+
+	call dgsder (gs, x, y, zfit, npts, 1, 1)
+	call printf ("\nZDER: 1 1\n")
+	do i = 1, npts {
+	    call printf ("%0.7g ")
+		call pargd (zfit[i])
+	    if (mod (i, 9) == 0)
+		call printf ("\n")
+	}
+	call printf ("\n")
+
+	# Refit the surface point by point.
+	call dgszero (gs)
+	do i = 1, npts {
+	    call dgsaccum (gs, x[i], y[i], z[i], w[i], WTS_UNIFORM)
+	}
+	if (stype > 3)
+	    call dgssolve1 (gs, ier)
+	else
+	    call dgssolve (gs, ier)
+	call printf ("\nACCUM FIT ERROR CODE: %d\n")
+	    call pargi (ier)
+	call dgsrej (gs, x[1], y[1], z[1], w[1], WTS_UNIFORM)
+	call dgsrej (gs, x[npts], y[npts], z[npts], w[npts], WTS_UNIFORM)
+	call dgsaccum (gs, x[1], y[1], z[1], w[1], WTS_UNIFORM)
+	call dgsaccum (gs, x[npts], y[npts], z[npts], w[npts], WTS_UNIFORM)
+	call dgssolve (gs, ier)
+	call printf ("\nREJ FIT ERROR CODE: %d\n")
+	    call pargi (ier)
+
+	call dgscoeff (gs, coeff, ncoeff)
+	call printf ("GSACCUM coeff:\n")
+	call printf ("first %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (gs, 1, 1))
+	    call pargd (dgsgcoeff (gs, xorder, 1))
+	do i = 1, ncoeff {
+	    call printf ("%d  %0.14g\n")
+		call pargi (i)
+		call pargd (coeff[i])
+	}
+	call printf ("last %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (gs, 1, yorder))
+	    call pargd (dgsgcoeff (gs, xorder, yorder))
+	call printf ("\n")
+
+	# Save and restore.
+	call dgssave (gs, save)
+	call dgsfree (gs)
+	call dgsrestore (gs, save)
+
+	call dgscoeff (gs, coeff, ncoeff)
+	call printf ("RESTORE coeff:\n")
+	call printf ("first %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (gs, 1, 1))
+	    call pargd (dgsgcoeff (gs, xorder, 1))
+	do i = 1, ncoeff {
+	    call printf ("%d  %0.14g\n")
+		call pargi (i)
+		call pargd (coeff[i])
+	}
+	call printf ("last %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (gs, 1, yorder))
+	    call pargd (dgsgcoeff (gs, xorder, yorder))
+	call printf ("\n")
+
+	# Add two surfaces.
+	call dgsadd (gs, gs, ags)
+	call dgscoeff (ags, coeff, ncoeff)
+	call printf ("GSADD coeff:\n")
+	call printf ("first %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (ags, 1, 1))
+	    call pargd (dgsgcoeff (ags, xorder, 1))
+	do i = 1, ncoeff {
+	    call printf ("%d  %0.14g\n")
+		call pargi (i)
+		call pargd (coeff[i])
+	}
+	call printf ("last %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (ags, 1, yorder))
+	    call pargd (dgsgcoeff (ags, xorder, yorder))
+	call printf ("\n")
+
+	# Subtract two surfaces.
+	call dgssub (gs, gs, sgs)
+	call dgscoeff (sgs, coeff, ncoeff)
+	call printf ("GSSUB coeff:\n")
+	call printf ("first %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (sgs, 1, 1))
+	    call pargd (dgsgcoeff (sgs, xorder, 1))
+	do i = 1, ncoeff {
+	    call printf ("%d  %0.14g\n")
+		call pargi (i)
+		call pargd (coeff[i])
+	}
+	call printf ("last %0.14g %0.14g\n")
+	    call pargd (dgsgcoeff (sgs, 1, yorder))
+	    call pargd (dgsgcoeff (sgs, xorder, yorder))
+	call printf ("\n")
+
+	call dgsfree (gs)
+	call dgsfree (ags)
+	call dgsfree (sgs)
+end