Initial commit

1 files changed, 986 insertions, 0 deletions

diff --git a/noao/onedspec/t_sfit.x b/noao/onedspec/t_sfit.x
new file mode 100644
index 00000000..54e896bb
--- /dev/null
+++ b/noao/onedspec/t_sfit.x
@@ -0,0 +1,986 @@
+include	<imhdr.h>
+include <pkg/gtools.h>
+include <pkg/rg.h>
+include <math/curfit.h>
+include	<error.h>
+include	<smw.h>
+
+# SFIT -- Fit a function to spectra and output the fit, difference or
+# ratio; or print the power series coefficients of the fit.  The fitting
+# parameters may be set interactively using the icfit package.
+
+# Image header keywords for saving the previous fit
+
+define	SFT_KW		"SFIT"
+define	SFT_KWB		"SFITB"
+
+# Choices for the type of output
+
+define	OUT_TYPES	"|data|fit|difference|ratio|"
+
+define	DATA		1
+define	FIT		2
+define	DIFFERENCE	3
+define	RATIO		4
+
+# Choices for the interactive prompts
+# (the 1st define is for clgwrd (strdic), the 2nd for CL enumeration)
+# Note that the CL assumes that the separator is `|'.
+
+define	SFT_ANS1	"|yes|no|skip|YES|NO|SKIP|"
+define	SFT_ANS1X	"yes|no|skip|YES|NO|SKIP"
+
+define	SFT_ANS2	"|spectrum|image|all|cancel|"
+define	SFT_ANS2X	"spectrum|image|all|cancel"
+
+define	LEN_ANS		7
+
+define	YES_ONCE	1
+define	NO_ONCE		2
+define	SKIP_ONCE	3
+define	YES_ALWAYS	4
+define	NO_ALWAYS	5
+define	SKIP_ALWAYS	6
+
+define	SKIP_SPEC	1
+define	SKIP_IMAGE	2
+define	SKIP_ALL	3
+define	SKIP_CANCEL	4
+
+# Switches and pointers
+
+define	SFT_OFF		22
+
+define	INTERACTIVE	Memi[$1]		# all spectra are noninteractive
+define	REPLACE		Memi[$1+1]		# replace rejected points?
+define	WAVESCALE	Memi[$1+2]		# X is wavelength if possible
+define	LOGSCALE	Memi[$1+3]		# axes are logarithmic
+define	OVERRIDE	Memi[$1+4]		# allow lines to be redone
+define	LISTONLY	Memi[$1+5]		# don't modify images
+define	OUTTYPE		Memi[$1+6]		# output type code
+
+define	GRAPH_OPEN	Memi[$1+7]		# keep track of gopen
+define	LOG_TO_STDOUT	Memi[$1+8]		# STDOUT/ERR is used
+define	PROMPT		Memi[$1+9]		# prompt flag
+define	QUIET		Memi[$1+10]		# quiet flag
+
+define	RGIN		Memi[$1+11]		# lines specified
+define	RGFIT		Memi[$1+12]		# all lines to fit
+define	RGREFIT		Memi[$1+13]		# those to fit again
+define	RGINB		Memi[$1+14]		# bands specified
+define	RGFITB		Memi[$1+15]		# all bands to fit
+define	RGREFITB	Memi[$1+16]		# those to fit again
+
+define	NLOGFD		Memi[$1+17]		# number of logfiles
+define	LOGFD		Memi[$1+18]		# array of logfiles
+
+define	IC		Memi[$1+19]		# current ic descriptor
+define	YMAX		Memi[$1+20]		# max number of lines
+define	BMAX		Memi[$1+21]		# max number of lines
+define	IC_DESC		Memi[$1+SFT_OFF+($3-1)*YMAX($1)+$2-1]	# ic descriptors
+
+
+# T_SFIT -- Entry point for the task.  Read parameters,
+# initialize structures and loop over the image templates.
+
+procedure t_sfit ()
+
+pointer	listin, listout, input, output, graphics
+pointer	sf, gp, gt, in, out, mw, sh, sp
+int	stat
+
+int	sft_icfit(), imtgetim(), gt_init(), imtlen()
+bool	clgetb()
+pointer	imtopenp()
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (graphics, SZ_FNAME, TY_CHAR)
+
+	# Open the image templates
+	listin = imtopenp ("input")
+
+	if (clgetb ("listonly"))
+	    listout = NULL
+	else {
+	    listout = imtopenp ("output")
+	    if (imtlen (listin) != imtlen (listout)) {
+		call imtclose (listin)
+		call imtclose (listout)
+		call sfree (sp)
+		call error (1, "Input and output image lists do not match")
+	    }
+	}
+
+	# Initialize the various descriptors
+	iferr (call sft_init (listin, listout, sf)) {
+	    call imtclose (listin)
+	    if (listout != NULL)
+		call imtclose (listout)
+	    call sfree (sp)
+	    call erract (EA_ERROR)
+	}
+
+	# The graphics pointers are passed explicitly
+	if (INTERACTIVE(sf) == YES) {
+	    call clgstr ("graphics", Memc[graphics], SZ_FNAME)
+	    gt = gt_init()
+	    call gt_sets (gt, GTTYPE, "line")
+	}
+
+	# Fit the lines in each input image.
+
+	while (imtgetim (listin, Memc[input], SZ_FNAME) != EOF) {
+
+	    if (listout != NULL)
+		stat = imtgetim (listout, Memc[output], SZ_FNAME)
+	    else
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+
+	    iferr {
+		call sft_immap (Memc[input], Memc[output],
+		    in, out, mw, sh, sf, gp)
+	    } then {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    stat = sft_icfit (in, out, mw, sh, sf, gp, gt, Memc[graphics])
+
+	    call sft_unmap (in, out, mw, sh, sf)
+
+	    if (stat == EOF)
+		break
+	}
+
+	if (INTERACTIVE(sf) == YES)
+	    call gt_free (gt)
+	if (GRAPH_OPEN(sf) == YES)
+	    call gclose (gp)
+
+	call sft_close (sf)
+	call imtclose (listin)
+	if (listout != NULL)
+	    call imtclose (listout)
+	call sfree (sp)
+end
+
+
+# SFT_INIT -- Initialize templates, ranges, logfiles, type, icfit descriptors.
+
+procedure sft_init (listin, listout, sf)
+
+pointer	listin, listout		#I Image template descriptors
+pointer	sf			#I Pointer to task switches
+
+pointer	input, output, im, sp, mw
+int	ymax, bmax, i, j
+
+real	clgetr()
+bool	clgetb()
+int	clgwrd(), clgeti(), btoi(), strlen()
+int	rg_next(), imtgetim(), imaccess(), xt_logopen()
+pointer	rg_ranges(), immap(), smw_openim()
+errchk	immap, smw_openim
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_LINE, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+
+	sf = NULL
+	iferr {
+
+	    # find the maximum number of lines and bands (spectra)
+	    ymax = 0
+	    bmax = 0
+
+	    while (imtgetim (listin, Memc[input], SZ_LINE) != EOF)
+		if (imaccess (Memc[input], READ_ONLY) == YES) {
+		    im = immap (Memc[input], READ_ONLY, 0)
+		    mw = smw_openim (im)
+		    ymax = max (SMW_LLEN(mw,2), ymax)
+		    bmax = max (SMW_LLEN(mw,3), bmax)
+		    call smw_close (mw)
+		    call imunmap (im)
+		}
+	    call imtrew (listin)
+
+	    if (listout != NULL) {
+		while (imtgetim (listout, Memc[output], SZ_FNAME) != EOF)
+		    if (imaccess (Memc[output], READ_ONLY) == YES) {
+			im = immap (Memc[output], READ_ONLY, 0)
+			mw = smw_openim (im)
+			ymax = max (SMW_LLEN(mw,2), ymax)
+			bmax = max (SMW_LLEN(mw,3), bmax)
+			call smw_close (mw)
+			call imunmap (im)
+		    }
+		call imtrew (listout)
+	    }
+
+	    # allocate space for the task switch structure
+	    call malloc (sf, SFT_OFF + ymax * bmax, TY_STRUCT)
+
+	    YMAX(sf) = ymax
+	    BMAX(sf) = bmax
+
+	    # NULL the pointers for error handling
+	    RGIN(sf)	= NULL
+	    RGINB(sf)	= NULL
+	    NLOGFD(sf)	= 0
+	    do j = 1, BMAX(sf)
+		do i = 1, YMAX(sf)
+		    IC_DESC(sf,i,j) = NULL
+
+	    # Set the switches
+	    INTERACTIVE(sf) = btoi (clgetb ("interactive"))
+	    REPLACE(sf) = btoi (clgetb ("replace"))
+	    WAVESCALE(sf) = btoi (clgetb ("wavescale"))
+	    LOGSCALE(sf) = btoi (clgetb ("logscale"))
+	    OVERRIDE(sf) = btoi (clgetb ("override"))
+	    LISTONLY(sf) = btoi (clgetb ("listonly"))
+	    GRAPH_OPEN(sf) = NO
+	    PROMPT(sf) = INTERACTIVE(sf)
+	    QUIET(sf) = btoi (INTERACTIVE(sf) == NO)
+
+	    # Expand the range specification, allow either hyphens or colons
+
+	    call clgstr ("lines", Memc[input], SZ_LINE)
+	    do i = 1, strlen (Memc[input])
+		if (Memc[input+i-1] == '-')
+		    Memc[input+i-1] =  ':'
+		else if (Memc[input+i-1] == 'x' || Memc[input+i-1] == 'X')
+		    call error (1, "Range step (`x' notation) not implemented")
+
+	    RGIN(sf) = rg_ranges (Memc[input], 1, YMAX(sf))
+	    call rg_order (RGIN(sf))
+	    call rg_merge (RGIN(sf))
+
+	    call clgstr ("bands", Memc[input], SZ_LINE)
+	    do i = 1, strlen (Memc[input])
+		if (Memc[input+i-1] == '-')
+		    Memc[input+i-1] =  ':'
+		else if (Memc[input+i-1] == 'x' || Memc[input+i-1] == 'X')
+		    call error (1, "Range step (`x' notation) not implemented")
+
+	    RGINB(sf) = rg_ranges (Memc[input], 1, BMAX(sf))
+	    call rg_order (RGINB(sf))
+	    call rg_merge (RGINB(sf))
+
+	    i = 0
+	    j = 0
+	    if (rg_next (RGIN(sf), i) == EOF || rg_next (RGINB(sf), j) == EOF)
+		call error (1, "With range specification for `lines or bands'")
+	    else {
+		# Open the initial icfit descriptor
+		call ic_open (IC(sf))
+
+		call clgstr ("sample", Memc[input], SZ_LINE)
+		call ic_pstr (IC(sf), "sample", Memc[input])
+		call clgstr ("function", Memc[input], SZ_LINE)
+		call ic_pstr (IC(sf), "function", Memc[input])
+
+		call ic_puti (IC(sf), "naverage", clgeti ("naverage"))
+		call ic_puti (IC(sf), "order", clgeti ("order"))
+		call ic_putr (IC(sf), "low", clgetr ("low_reject"))
+		call ic_putr (IC(sf), "high", clgetr ("high_reject"))
+		call ic_puti (IC(sf), "niterate", clgeti ("niterate"))
+		call ic_putr (IC(sf), "grow", clgetr ("grow"))
+		call ic_puti (IC(sf), "markrej", btoi (clgetb ("markrej")))
+
+		IC_DESC(sf,i,j) = IC(sf)
+	    }
+
+	    # Get the desired output type
+	    OUTTYPE(sf) = clgwrd ("type", Memc[input], SZ_LINE, OUT_TYPES)
+
+	    # Open the logfiles
+	    NLOGFD(sf) = xt_logopen ("logfiles", "SFIT:", LOGFD(sf),
+		LOG_TO_STDOUT(sf))
+
+	} then {
+	    call sfree (sp)
+	    call sft_close (sf)
+	    call erract (EA_ERROR)
+	}
+
+	call sfree (sp)
+	return
+end
+
+
+# SFT_CLOSE -- Close the various descriptors.
+
+procedure sft_close (sf)
+
+pointer	sf			#I Pointer to task switches
+
+int	i, j
+
+begin
+	if (sf != NULL) {
+	    if (RGIN(sf) != NULL)
+		call rg_free (RGIN(sf))
+	    if (RGINB(sf) != NULL)
+		call rg_free (RGINB(sf))
+	    if (NLOGFD(sf) != 0)
+		call xt_logclose (LOGFD(sf), NLOGFD(sf), "END:")
+	    do j = 1, BMAX(sf)
+		do i = 1, YMAX(sf)
+		    if (IC_DESC(sf,i,j) != NULL)
+			call ic_closer (IC_DESC(sf,i,j))
+	    call mfree (sf, TY_STRUCT)
+	}
+end
+
+
+# SFT_IMMAP -- Map images for sfit.
+
+procedure sft_immap (input, output, in, out, mw, sh, sf, gp)
+
+char	input[ARB]		#I Input image name
+char	output[ARB]		#I Output image name
+pointer	in, out			#O IMIO pointers
+pointer	mw			#O MWCS pointer
+pointer	sh			#O SHDR pointer
+pointer	sf			#I Pointer for task switches
+pointer	gp			#I GIO pointer
+
+int	i, ax1, ax2, ax3
+pointer	inroot, insect, outroot, outsect, b1, b2
+pointer	sp, inranges, outranges
+pointer	rgin, rgout, rgtmp, rgtmpb
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+char	emsg[SZ_LINE]
+
+int	imaccess(), imaccf(), imgnlr(), impnlr(), strcmp()
+pointer	immap(), smw_openim()
+pointer	rg_ranges(), rg_window(), rg_union(), rg_intersect()
+errchk	immap, smw_openim
+
+define	err_	13
+
+begin
+	call smark (sp)
+	call salloc (inroot, SZ_FNAME, TY_CHAR)
+	call salloc (insect, SZ_FNAME, TY_CHAR)
+	call salloc (outroot, SZ_FNAME, TY_CHAR)
+	call salloc (outsect, SZ_FNAME, TY_CHAR)
+	call salloc (inranges, SZ_LINE, TY_CHAR)
+	call salloc (outranges, SZ_LINE, TY_CHAR)
+
+	in = NULL
+	out = NULL
+	mw = NULL
+	sh = NULL
+	RGFIT(sf) = NULL
+	RGREFIT(sf) = NULL
+	RGFITB(sf) = NULL
+	RGREFITB(sf) = NULL
+
+	call imgimage (input, Memc[inroot], SZ_FNAME)
+	call imgsection (input, Memc[insect], SZ_FNAME)
+
+	call imgimage (output, Memc[outroot], SZ_FNAME)
+	call imgsection (output, Memc[outsect], SZ_FNAME)
+
+	if (Memc[insect] != EOS || Memc[outsect] != EOS) {
+
+	    call sprintf (emsg, SZ_LINE, "Sections not allowed (%s --> %s)")
+		call pargstr (input)
+		call pargstr (output)
+	    goto err_
+
+	} else if (imaccess (Memc[inroot], READ_ONLY) == NO) {
+
+	    call sprintf (emsg, SZ_LINE, "Cannot access %s")
+		call pargstr (input)
+	    goto err_
+
+	} else if (LISTONLY(sf) == YES) {
+
+	    # The `out = in' allows the ranges code at the end of this
+	    # procedure to cover all cases (with a little inefficiency).
+	    # No check on the sizes of the input and output images.
+
+	    in = immap (Memc[inroot], READ_ONLY, 0)
+	    out = in
+
+	} else if (strcmp (Memc[inroot], Memc[outroot]) == 0) {
+
+	    # Overwrite the input image.
+	    in = immap (Memc[inroot], READ_WRITE, 0)
+	    out = in
+
+	} else if (imaccess (Memc[outroot], READ_WRITE) == NO) {
+
+	    in = immap (Memc[inroot], READ_ONLY, 0)
+	    out = immap (Memc[outroot], NEW_COPY, in)
+	    if (IM_PIXTYPE(out) != TY_DOUBLE)
+		IM_PIXTYPE(out) = TY_REAL
+
+	    # Do this since imcopy is unimplemented
+
+	    call amovkl (long(1), v1, IM_MAXDIM)
+	    call amovkl (long(1), v2, IM_MAXDIM)
+
+	    while (imgnlr (in, b1, v1) != EOF && impnlr (out, b2, v2) != EOF)
+		call amovr (Memr[b1], Memr[b2], IM_LEN(in, 1))
+
+	} else {
+
+	    in = immap (Memc[inroot], READ_ONLY, 0)
+	    out = immap (Memc[outroot], READ_WRITE, 0)
+
+	    # This relies on the axes beyond IM_NDIM(im) being unity
+
+	    do i = 1, max (IM_NDIM(in), IM_NDIM(out))
+		if (IM_LEN(in, i) != IM_LEN(out, i)) {
+		    call sprintf (emsg, SZ_LINE, "%s & %s aren't the same size")
+			call pargstr (Memc[inroot])
+			call pargstr (Memc[outroot])
+		    goto err_
+		}
+
+	}
+
+	do i = 4, IM_NDIM(in)
+	    if (IM_LEN(in, i) != 1) {
+		call sprintf (emsg, SZ_LINE, "Too many dimensions for %s")
+		    call pargstr (Memc[inroot])
+		goto err_
+	    }
+
+	if (imaccf (in, SFT_KW) == YES)
+	    call imgstr (in, SFT_KW, Memc[inranges], SZ_LINE)
+	else
+	    call strcpy ("", Memc[inranges], SZ_LINE)
+
+	if (imaccf (out, SFT_KW) == YES)
+	    call imgstr (out, SFT_KW, Memc[outranges], SZ_LINE)
+	else {
+	    call strcpy ("", Memc[outranges], SZ_LINE)
+	    call imastr (out, SFT_KW, Memc[outranges])
+	}
+
+	mw = smw_openim (in)
+	ax1 = SMW_LLEN(mw,1)
+	ax2 = SMW_LLEN(mw,2)
+	ax3 = SMW_LLEN(mw,3)
+
+	rgin = rg_ranges (Memc[inranges], 1, ax2)
+	rgout = rg_ranges (Memc[outranges], 1, ax2)
+	rgtmp = rg_union (rgin, rgout)
+	call rg_free (rgin)
+	call rg_free (rgout)
+
+	if (imaccf (in, SFT_KWB) == YES)
+	    call imgstr (in, SFT_KWB, Memc[inranges], SZ_LINE)
+	else
+	    call strcpy ("", Memc[inranges], SZ_LINE)
+
+	if (imaccf (out, SFT_KWB) == YES)
+	    call imgstr (out, SFT_KWB, Memc[outranges], SZ_LINE)
+	else {
+	    call strcpy ("", Memc[outranges], SZ_LINE)
+	    call imastr (out, SFT_KWB, Memc[outranges])
+	}
+
+	rgin = rg_ranges (Memc[inranges], 1, ax3)
+	rgout = rg_ranges (Memc[outranges], 1, ax3)
+	rgtmpb = rg_union (rgin, rgout)
+	call rg_free (rgin)
+	call rg_free (rgout)
+
+	if (OVERRIDE(sf) == YES) {
+	    RGFIT(sf) = rg_window (RGIN(sf), 1, ax2)
+	    RGREFIT(sf) = rgtmp
+	    RGFITB(sf) = rg_window (RGINB(sf), 1, ax3)
+	    RGREFITB(sf) = rgtmpb
+	} else {
+	    call rg_inverse (rgtmp, 1, ax2)
+	    RGFIT(sf) = rg_intersect (RGIN(sf), rgtmp)
+	    RGREFIT(sf) = rg_ranges ("0", 1, 2)
+	    call rg_free (rgtmp)
+	    #call rg_inverse (rgtmpb, 1, ax3)
+	    #RGFITB(sf) = rg_intersect (RGINB(sf), rgtmpb)
+	    #RGREFITB(sf) = rg_ranges ("0", 1, 2)
+	    #call rg_free (rgtmpb)
+	    RGFITB(sf) = rg_window (RGINB(sf), 1, ax3)
+	    RGREFITB(sf) = rgtmpb
+	}
+
+	if (RG_NPTS(RGFIT(sf)) <= 0) {
+	    call sprintf (emsg, SZ_LINE, "No lines left to fit for %s")
+		call pargstr (Memc[inroot])
+	    goto err_
+	}
+
+	call sfree (sp)
+	return
+
+err_	call sfree (sp)
+	call sft_unmap (in, out, mw, sh, sf)
+	if (GRAPH_OPEN(sf) == YES) {
+	    call gclose (gp)
+	    GRAPH_OPEN(sf) = NO
+	}
+	# STDERR should get flushed AFTER closing graphics
+	call error (1, emsg)
+end
+
+
+# SFT_UNMAP -- Unmap images for sfit.
+
+procedure sft_unmap (in, out, mw, sh, sf)
+
+pointer	in, out			#I IMIO pointers
+pointer	mw			#I MWCS pointer
+pointer	sh			#I SHDR pointer
+pointer	sf			#I Task structure pointer
+
+begin
+	call shdr_close (sh)
+	if (mw != NULL)
+	    call smw_close (mw)
+	if (out != NULL && out != in)
+	    call imunmap (out)
+	if (in != NULL)
+	    call imunmap (in)
+	if (RGFIT(sf) != NULL)
+	    call rg_free (RGFIT(sf))
+	if (RGREFIT(sf) != NULL)
+	    call rg_free (RGREFIT(sf))
+	if (RGFITB(sf) != NULL)
+	    call rg_free (RGFITB(sf))
+	if (RGREFITB(sf) != NULL)
+	    call rg_free (RGREFITB(sf))
+end
+
+
+# SFT_ICFIT -- Given the image descriptors determine the fitting function
+# for each line and output the fit, difference, ratio or coefficients.
+
+int procedure sft_icfit (in, out, mw, sh, sf, gp, gt, graphics)
+
+pointer	in, out				#I IMIO pointers
+pointer	mw				#I MWCS pointer
+pointer	sh				#I SHDR pointer
+pointer	sf				#I Pointer for task switches
+pointer	gp				#I GIO pointer
+pointer	gt				#I GTOOLS pointer
+char	graphics[ARB]			#I Graphics device
+
+pointer	sp, wts, cv, data
+int	line, band, i, j, n
+
+int	sft_getline()
+pointer	gopen(), imps3r()
+real	sft_efncr()
+extern	sft_efncr
+
+begin
+	call smark (sp)
+	call salloc (wts, SMW_LLEN(mw,1), TY_REAL)
+
+	line = 0
+	band = 0
+	while (sft_getline (in, mw, sh, sf, gt, line, band) != EOF) {
+
+	    call amovkr (1., Memr[wts], SN(sh))
+
+	    if (QUIET(sf) == NO) {
+		if (GRAPH_OPEN(sf) == NO) {
+		    gp = gopen (graphics, NEW_FILE, STDGRAPH)
+		    GRAPH_OPEN(sf) = YES
+		}
+		call icg_fit (IC(sf), gp, "cursor", gt, cv, Memr[SX(sh)],
+		    Memr[SY(sh)], Memr[wts], SN(sh))
+	    } else
+		call ic_fit (IC(sf), cv, Memr[SX(sh)], Memr[SY(sh)],
+		    Memr[wts], SN(sh), YES, YES, YES, YES)
+
+	    if (LISTONLY(sf) == NO) {
+		i = LINDEX(sh,1)
+		j = LINDEX(sh,2)
+		n = SMW_LLEN(mw,1)
+		switch (SMW_LAXIS(mw,1)) {
+		case 1:
+		    data = imps3r (out, 1, n, i, i, j, j)
+		case 2:
+		    data = imps3r (out, i, i, 1, n, j, j)
+		case 3:
+		    data = imps3r (out, i, i, j, j,  1, n)
+		}
+		if (SN(sh) < n)
+		    call aclrr (Memr[data], n)
+
+		switch (OUTTYPE(sf)) {
+		case DATA:
+		    if (REPLACE(sf) == YES)
+			call ic_clean (IC(sf), cv, Memr[SX(sh)], Memr[SY(sh)],
+			    Memr[wts], SN(sh))
+		    call amovr (Memr[SY(sh)], Memr[data], SN(sh))
+		    call sft_update (out, mw, line, band)
+		case FIT:
+		    call cvvector (cv, Memr[SX(sh)], Memr[data], SN(sh))
+		    call sft_update (out, mw, line, band)
+		case DIFFERENCE:
+		    call cvvector (cv, Memr[SX(sh)], Memr[data], SN(sh))
+		    if (REPLACE(sf) == YES)
+			call ic_clean (IC(sf), cv, Memr[SX(sh)], Memr[SY(sh)],
+			    Memr[wts], SN(sh))
+		    call asubr (Memr[SY(sh)], Memr[data], Memr[data], SN(sh))
+		    call sft_update (out, mw, line, band)
+		case RATIO:
+		    call cvvector (cv, Memr[SX(sh)], Memr[data], SN(sh))
+		    if (REPLACE(sf) == YES)
+			call ic_clean (IC(sf), cv, Memr[SX(sh)], Memr[SY(sh)],
+			    Memr[wts], SN(sh))
+		    call advzr (Memr[SY(sh)], Memr[data], Memr[data], SN(sh),
+			sft_efncr)
+		    call sft_update (out, mw, line, band)
+		default:
+		    call error (1, "bad switch in sft_icfit")
+		}
+	    }
+
+	    call sft_power (in, line, cv, gp, sf)
+	    call cvfree (cv)
+
+	}
+
+	# This terminates the cursor (GIN) mode echoplex suppression in
+	# case the next sft_immap generates a password prompt from ZFIOKS.
+	# Note that any such password prompt (from the kernel!) will
+	# now show up on the status line, not the graphics plane.
+
+	if (GRAPH_OPEN(sf) == YES) {
+	    call printf ("\r")
+	    call flush (STDOUT)
+	}
+
+	call sfree (sp)
+	return (line)
+end
+
+
+# SFT_GETLINE -- Get image data to be fit.  Returns the line and band numbers.
+# Returns EOF when done.
+
+int procedure sft_getline (in, mw, sh, sf, gt, line, band)
+
+pointer	in			#I IMIO pointer
+pointer	mw			#I MWCS pointer
+pointer	sh			#I SHDR pointer
+pointer	sf			#I Pointer for task switches
+pointer	gt			#I GTOOLS pointer
+int	line			#U Line number
+int	band			#U Band number
+
+int	i
+bool	waveok
+char	ask[LEN_ANS]
+pointer	linebuf, rg1, rg2, sp
+
+int	clgwrd(), rg_next(), rg_inrange()
+pointer	rg_ranges(), rg_intersect()
+real	sft_efncr()
+extern	sft_efncr
+errchk	shdr_open
+
+define	again_		99
+
+begin
+	call smark (sp)
+	call salloc (linebuf, SZ_LINE, TY_CHAR)
+
+	if (band == 0)
+	    if (rg_next (RGFITB(sf), band) == EOF)
+		return (EOF)
+
+again_	if (rg_next (RGFIT(sf), line) == EOF) {
+	    line = 0
+	    if (rg_next (RGFITB(sf), band) == EOF)
+		return (EOF)
+	    goto again_
+	}
+
+	if (PROMPT(sf) == YES) {
+	    call clprintf ("ask.p_min", "%s")
+		call pargstr (SFT_ANS1X)
+
+	    if (rg_inrange (RGREFIT(sf), line) == YES &&
+		rg_inrange (RGREFITB(sf), band) == YES) {
+		call clprintf ("ask.p_prompt",
+		    "Refit [%d,%d] of %s w/ graph?        ")
+	    } else {
+		call clprintf ("ask.p_prompt",
+		    "Fit [%d,%d] of %s w/ graph?          ")
+	    }
+	    call pargi (line)
+	    call pargi (band)
+	    call pargstr (IM_HDRFILE(in))
+
+	    switch (clgwrd ("ask", ask, LEN_ANS, SFT_ANS1)) {
+
+	    case YES_ONCE:
+		QUIET(sf) = NO
+
+	    case NO_ONCE:
+		QUIET(sf) = YES
+
+	    case SKIP_ONCE:
+		goto again_
+
+	    case YES_ALWAYS:
+		QUIET(sf) = NO
+		PROMPT(sf) = NO
+
+	    case NO_ALWAYS:
+		QUIET(sf) = YES
+		PROMPT(sf) = NO
+
+	    case SKIP_ALWAYS:
+		call clprintf ("ask", "cancel")
+		call clprintf ("ask.p_min", "%s")
+		    call pargstr (SFT_ANS2X)
+		call clprintf ("ask.p_prompt",
+		    "Skip what?    (`all' to exit task)   ")
+
+		switch (clgwrd ("ask", ask, LEN_ANS, SFT_ANS2)) {
+
+		case SKIP_SPEC:
+		    call clprintf ("ask", "yes")
+		    # delete the spectrum from the list
+		    call sprintf (Memc[linebuf], SZ_LINE, "%d")
+			call pargi (line)
+
+		    rg1 = rg_ranges (Memc[linebuf], 1, SMW_LLEN(mw,2))
+		    call rg_inverse (rg1, 1, SMW_LLEN(mw,2))
+		    rg2 = rg_intersect (RGIN(sf), rg1)
+		    call rg_free (rg1)
+		    call rg_free (RGIN(sf))
+
+		    RGIN(sf) = rg2
+		    goto again_
+
+		case SKIP_IMAGE:
+		    call clprintf ("ask", "yes")
+		    return (EOF)
+
+		case SKIP_ALL:
+		    call clprintf ("ask", "yes")
+		    return (EOF)
+
+		case SKIP_CANCEL:
+		    call clprintf ("ask", "yes")
+		    line = line - 1
+		    goto again_
+
+		default:
+		    call error (1, "bad switch (2) in sft_getline")
+
+		}
+
+	    default:
+		call error (1, "bad switch (1) in sft_getline")
+
+	    }
+
+	}
+
+	call shdr_open (in, mw, line, band, INDEFI, SHDATA, sh)
+
+	if (LOGSCALE(sf) == YES)
+	    call alogr (Memr[SY(sh)], Memr[SY(sh)], SN(sh), sft_efncr)
+
+	if (WAVESCALE(sf) == YES) {
+	    waveok = true
+	} else
+	    waveok = false
+
+	if (!waveok)
+	    do i = 1, SN(sh)
+		Memr[SX(sh)+i-1] = i
+
+	if (LOGSCALE(sf) == YES)
+	    call alogr (Memr[SX(sh)], Memr[SX(sh)], SN(sh), sft_efncr)
+
+	# Initialize and/or update the icfit descriptor
+
+	if (IC_DESC(sf,line,band) == NULL) {
+	    call ic_open (IC_DESC(sf,line,band))
+	    call ic_copy (IC(sf), IC_DESC(sf,line,band))
+	    #call ic_pstr (IC_DESC(sf,line,band), "sample", "*")
+	}
+
+	IC(sf) = IC_DESC(sf,line,band)
+
+	call ic_putr (IC(sf), "xmin", min (Memr[SX(sh)], Memr[SX(sh)+SN(sh)-1]))
+	call ic_putr (IC(sf), "xmax", max (Memr[SX(sh)], Memr[SX(sh)+SN(sh)-1]))
+
+	if (QUIET(sf) == NO) {
+	    if (waveok && LOGSCALE(sf) == YES) {
+		call ic_pstr (IC(sf), "xlabel", "log wavelength")
+		call ic_pstr (IC(sf), "ylabel", "log data")
+	    } else if (LOGSCALE(sf) == YES) {
+		call ic_pstr (IC(sf), "xlabel", "log column")
+		call ic_pstr (IC(sf), "ylabel", "log data")
+	    } else if (waveok) {
+		call ic_pstr (IC(sf), "xlabel", "wavelength")
+		call ic_pstr (IC(sf), "ylabel", "")
+	    } else {
+		call ic_pstr (IC(sf), "xlabel", "column")
+		call ic_pstr (IC(sf), "ylabel", "")
+	    }
+
+	    call sprintf (Memc[linebuf], SZ_LINE, "%s, [%d,%d]\n%s")
+		call pargstr (IM_HDRFILE(in))
+		call pargi (line)
+		call pargi (band)
+		call pargstr (TITLE(sh))
+
+	    call gt_sets (gt, GTTITLE, Memc[linebuf])
+	}
+
+	call sfree (sp)
+	return (OK)
+end
+
+
+# SFT_EFNCR -- Called by advzr on division by zero or by alogr for a
+# zero or negative argument.
+
+real procedure sft_efncr (x)
+
+real	x
+
+begin
+	return (0.)
+end
+
+
+# SFT_POWER -- Transform the curfit output into a power series and
+# print the coefficients to the logfiles.  This should be modified to
+# print the errors as well.  That requires modifying the curfit routine
+# cvpower to deal with errors; and adding an icfit routine (or include
+# file define) that allows access to the dynamic arrays of sample points
+# that are initialized if the sample is less than the whole set of points.
+
+procedure sft_power (im, line, cv, gp, sf)
+
+pointer	im			#I IMIO descriptor for labeling
+int	line			#I Image line number for labeling
+pointer	cv			#I CURFIT pointer
+pointer	gp			#I GIO pointer for tidy output
+pointer sf			#I Pointer for task switches
+
+pointer	ps_coeff, linebuf, sp
+int	ncoeffs, i, j, fd
+
+int	cvstati(), strcmp()
+
+begin
+	if (NLOGFD(sf) <= 0)
+	    return
+
+	call smark (sp)
+	call salloc (linebuf, SZ_LINE, TY_CHAR)
+
+	# cvpower only works with legendre or chebyshev functions
+
+	call ic_gstr (IC(sf), "function", Memc[linebuf], SZ_LINE)
+	if (strcmp (Memc[linebuf], "legendre") != 0 &&
+	    strcmp (Memc[linebuf], "chebyshev") != 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	if (GRAPH_OPEN(sf) == YES && LOG_TO_STDOUT(sf) == YES) {
+	    call gclose (gp)
+	    GRAPH_OPEN(sf) = NO
+	}
+
+	ncoeffs = cvstati (cv, CVNCOEFF)
+	call salloc (ps_coeff, ncoeffs, TY_REAL)
+	call cvpower (cv, Memr[ps_coeff], ncoeffs)
+
+	do i = 1, NLOGFD(sf) {
+	    fd = Memi[LOGFD(sf)+i-1]
+
+	    call fprintf (fd, "Line %d of %s:\n\n")
+		call pargi (line)
+		call pargstr (IM_HDRFILE(im))
+
+	    call fprintf (fd, "      coeff           value\n")
+
+	    do j = 1, ncoeffs {
+		call fprintf (fd, "\t%d\t%12.5e\n")
+		    call pargi (j)
+		    call pargr (Memr[ps_coeff+j-1])
+	    }
+
+	    call fprintf (fd, "\n")
+	    call flush (fd)
+	}
+
+	call sfree (sp)
+end
+
+
+# SFT_UPDATE -- Update the keyword with completed spectrum.  Flush the pixels.
+
+procedure sft_update (im, mw, line, band)
+
+pointer	im			#I IMIO pointer
+pointer	mw			#I MWCS pointer
+int	line			#I Line just completed
+int	band			#I Band just completed
+
+pointer linebuf, rg1, rg2, rgold, sp
+
+pointer	rg_ranges(), rg_union()
+
+begin
+	call smark (sp)
+	call salloc (linebuf, SZ_LINE, TY_CHAR)
+
+	# this could be recoded to use "rg_add"
+
+	call sprintf (Memc[linebuf], SZ_LINE, "%d")
+	    call pargi (line)
+	rg1 = rg_ranges (Memc[linebuf], 1, SMW_LLEN(mw,2))
+
+	call imgstr (im, SFT_KW, Memc[linebuf], SZ_LINE)
+	rg2 = rg_ranges (Memc[linebuf], 1, SMW_LLEN(mw,2))
+
+	rgold = rg_union (rg1, rg2)
+	call rg_encode (rgold, Memc[linebuf], SZ_LINE)
+	call impstr (im, SFT_KW, Memc[linebuf])
+
+	call rg_free (rg1)
+	call rg_free (rg2)
+	call rg_free (rgold)
+
+	call sprintf (Memc[linebuf], SZ_LINE, "%d")
+	    call pargi (band)
+	rg1 = rg_ranges (Memc[linebuf], 1, SMW_LLEN(mw,3))
+
+	call imgstr (im, SFT_KWB, Memc[linebuf], SZ_LINE)
+	rg2 = rg_ranges (Memc[linebuf], 1, SMW_LLEN(mw,3))
+
+	rgold = rg_union (rg1, rg2)
+	call rg_encode (rgold, Memc[linebuf], SZ_LINE)
+	call impstr (im, SFT_KWB, Memc[linebuf])
+
+	call rg_free (rg1)
+	call rg_free (rg2)
+	call rg_free (rgold)
+
+	call imflush (im)
+	call sfree (sp)
+end