Initial commit

1 files changed, 232 insertions, 0 deletions

diff --git a/noao/onedspec/t_sinterp.x b/noao/onedspec/t_sinterp.x
new file mode 100644
index 00000000..2275a42b
--- /dev/null
+++ b/noao/onedspec/t_sinterp.x
@@ -0,0 +1,232 @@
+include	<imhdr.h>
+include	<math/curfit.h>
+
+# Interpolation mode
+define	SI_LINEAR		1
+define	SI_CURVES		2
+define	SI_LEGENDRE		3
+define	SI_CHEBYSHEV		4
+define	SI_SPLINE3		5
+define	SI_SPLINE1		6
+
+# T_SINTERP -- Interpolate for values in a table and optionally generate
+#		a spectral image
+#
+# A table of x,y pairs contained in a file is used to
+# find interpolated values, y, for any other given independent
+# variable, x. Extrapolation is performed if necessary.
+#
+# A series of values may be generated to generate a fine grid
+# through a coarse sampling for purposes of plotting. This is
+# done by setting the hidden parameter curve_gen to yes.
+# The starting point, ending point, and sampling interval
+# are also needed in this case (x1, x2, dx).
+#
+# If only a small number of values are needed to be interpolated
+# from the table, the user may enter a number of x's from either
+# a file or STDIN.
+
+procedure t_sinterp()
+
+real	x, y, x1, x2, dx
+int	npts, i
+int	filelist, tbl, in
+int	user_mode, imlen, order, maxpts
+char	fname[SZ_FNAME], tbl_file[SZ_FNAME]
+char	image[SZ_FNAME]
+char	interp[SZ_LINE]
+bool	gen, make_image
+pointer	im, pix, sp, xtab, ytab, cv
+
+int	clpopni(), clgfil(), open(), fscan(), nscan()
+int	clgeti(), clgwrd()
+real	clgetr()
+bool	clgetb()
+pointer	immap(), impl1r()
+
+begin
+	# Initialize interpolator
+	call intrp0 (1)
+	cv = NULL
+
+	# File containing x,y pairs in a table
+	call clgstr ("tbl_file", tbl_file, SZ_FNAME)
+
+	# Open table file and read as many points as possible
+	tbl = open (tbl_file, READ_ONLY, TEXT_FILE)
+
+	npts = 0
+	maxpts = clgeti ("tbl_size")
+
+	call smark (sp)
+	call salloc (xtab, maxpts, TY_REAL)
+	call salloc (ytab, maxpts, TY_REAL)
+
+	while (fscan(tbl) != EOF) {
+	    npts = npts + 1
+	    if (npts > maxpts)
+		call error (1, "Maximum table size exceeded.")
+	    call gargr (Memr[xtab+npts-1])
+	    call gargr (Memr[ytab+npts-1])
+	    if (nscan() < 2) {
+#		call eprintf ("Error reading x,y pairs\n")
+		npts = npts - 1
+	    }
+	}
+
+	call close (tbl)
+
+	if (npts < 1)
+	    call error (1, "Table has no entries.")
+
+	# Linear, spline, or CURFIT option interpolator?
+	user_mode = clgwrd ("interp_mode", interp, SZ_LINE,
+	    ",linear,curves,legendre,chebyshev,spline3,spline1")
+
+	if (user_mode > 2 && user_mode <= 6)
+	    order = clgeti ("order")
+
+	# Generate a curve?
+	gen = clgetb ("curve_gen")
+
+	# Or an image?
+	make_image = clgetb ("make_image")
+
+	if (gen || make_image) {
+	    x1 = clgetr ("x1")
+	    x2 = clgetr ("x2")
+	    dx = clgetr ("dx")
+	    imlen = clgeti ("npts")
+
+	    # The above four variables overdefine the function
+	    # One (other than x1) must be 0.0 --> solve for it
+	    if (x2 == 0.0)
+		x2 = x1 + (imlen-1) * dx
+	    else if (dx == 0.0)
+		dx = (x2 - x1) / (imlen - 1)
+
+	    imlen = nint ((x2 - x1) / dx + 1)
+
+	    # Verify that dx will not cause an infinite loop
+	    if (dx == 0.0 || dx * (x2-x1) < 0.0)
+		call error (1, "Interval paramater dx implies infinite loop.")
+
+	    if (make_image) {
+		call clgstr ("image", image, SZ_FNAME)
+		im = immap (image, NEW_IMAGE, 0)
+
+		IM_NDIM (im) = 1
+		IM_LEN (im, 1) = imlen
+		IM_PIXTYPE (im) = TY_REAL
+
+		pix = impl1r (im)
+
+		do i = 1, imlen {
+		    x = x1 + (i - 1) * dx
+		    call gen_pixel (Memr[xtab], Memr[ytab], npts, 
+			user_mode, order, x, y, cv)
+		    Memr[pix+i-1] = y
+		}
+
+		call imaddr (im, "CRVAL1", x1)
+		call imaddr (im, "CDELT1", dx)
+		call imaddr (im, "CD1_1", dx)
+		call imaddr (im, "CRPIX1", 1.)
+		call imaddi (im, "DC-FLAG", 0)
+
+		call imunmap (im)
+	    } else {
+		do i = 1, imlen {
+		    x = x1 + (i - 1) * dx
+		    call gen_pixel (Memr[xtab], Memr[ytab], npts, 
+			user_mode, order, x, y, cv)
+		    call printf ("%12.5g  %12.5g\n")
+		        call pargr (x)
+		        call pargr (y)
+		}
+		call flush (STDOUT)
+	    }
+
+	# No, just one point at a time
+	} else {
+
+	    # Open input list
+	    filelist = clpopni ("input")
+
+	    while (clgfil (filelist, fname, SZ_FNAME) != EOF) {
+		in = open (fname, READ_ONLY, TEXT_FILE)
+
+		# Process input requests
+		while (fscan(in) != EOF) {
+		    call gargr (x)
+
+		    call gen_pixel (Memr[xtab], Memr[ytab], npts, 
+			user_mode, order, x, y, cv)
+		    call printf ("%12.5g  %12.5g\n")
+			call pargr (x)
+			call pargr (y)
+		    call flush (STDOUT)
+		}
+
+		call close (in)
+	    }
+
+	    call clpcls (filelist)
+	}
+
+	call cvfree (cv)
+	call sfree (sp)
+end
+
+# GEN_PIXEL -- Generate a pixel value using specified interpolator
+
+procedure gen_pixel (xtab, ytab, npts, mode, order, x, y, cv)
+
+real	xtab[ARB], ytab[ARB]
+int	npts
+real	x
+int	mode, order
+real	y
+pointer	cv
+
+int	fit, ier
+pointer	wt, sp
+
+real	cveval()
+
+begin
+	# Interpolate after selecting option
+	switch (mode) {
+	    case SI_LINEAR:
+		call lintrp (1, xtab, ytab, npts, x, y, ier)
+
+	    case SI_CURVES:
+		call intrp  (1, xtab, ytab, npts, x, y, ier)
+
+	    default:
+		if (cv == NULL) {
+		    call smark (sp)
+		    call salloc (wt, npts, TY_REAL)
+		    call amovkr (1.0, Memr[wt], npts)
+
+		    switch (mode) {
+	    		case SI_LEGENDRE:
+			    fit = LEGENDRE
+	    		case SI_CHEBYSHEV:
+			    fit = CHEBYSHEV
+	    		case SI_SPLINE3:
+			    fit = SPLINE3
+	    		case SI_SPLINE1:
+			    fit = SPLINE1
+			default:
+			    fit = SPLINE1
+		    }
+
+		    call cvinit (cv, fit, order, xtab[1], xtab[npts])
+		    call cvfit (cv, xtab, ytab, Memr[wt], npts, WTS_UNIFORM,
+			ier)
+		    call sfree (sp)
+		}
+		y = cveval (cv, x)
+	}
+end