1 files changed, 310 insertions, 0 deletions
diff --git a/noao/imred/ccdred/src/setoverscan.x b/noao/imred/ccdred/src/setoverscan.x
new file mode 100644
index 00000000..e344aa92
--- /dev/null
+++ b/noao/imred/ccdred/src/setoverscan.x
@@ -0,0 +1,310 @@
+include	<imhdr.h>
+include	<imset.h>
+include <pkg/gtools.h>
+include	<pkg/xtanswer.h>
+include	"ccdred.h"
+
+
+# SET_OVERSCAN -- Set the overscan vector.
+#
+#   1.  Return immediately if the overscan correction is not requested or
+#	if the image has been previously corrected.
+#   2.	Determine the overscan columns or lines.  This may be specifed
+#	directly or indirectly through the image header or symbol table.
+#   3.  Determine the type of overscan.
+#   4.	If fitting the overscan average the overscan columns or lines and
+#     	fit a function with the ICFIT routines to smooth the overscan vector.
+#   5.  Set the processing flag.
+#   6.  Log the operation (to user, logfile, and output image header).
+
+procedure set_overscan (ccd)
+
+pointer	ccd			# CCD structure pointer
+
+int	i, first, last, navg, npts, type
+int	nc, nl, c1, c2, l1, l2
+pointer	sp, str, errstr, func, buf, x, overscan
+
+int	clgwrd()
+real	asumr()
+bool	clgetb(), ccdflag()
+pointer	imgl2r(), imgs2r()
+errchk	imgl2r, imgs2r, fit_overscan
+
+begin
+	# Check if the user wants this operation or if it has been done.
+	if (!clgetb ("overscan") || ccdflag (IN_IM(ccd), "overscan"))
+	    return
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (errstr, SZ_LINE, TY_CHAR)
+	call salloc (func, SZ_LINE, TY_CHAR)
+	call imstats (IN_IM(ccd), IM_IMAGENAME, Memc[str], SZ_LINE)
+
+	# Check bias section.
+	nc = IM_LEN(IN_IM(ccd),1)
+	nl = IM_LEN(IN_IM(ccd),2)
+	c1 = BIAS_C1(ccd)
+	c2 = BIAS_C2(ccd)
+	l1 = BIAS_L1(ccd)
+	l2 = BIAS_L2(ccd)
+	if ((c1 < 1) || (c2 > nc) || (l1 < 1) || (l2 > nl)) {
+	    call sprintf (Memc[errstr], SZ_LINE,
+		"Error in bias section: image=%s[%d,%d], biassec=[%d:%d,%d:%d]")
+		call pargstr (Memc[str])
+		call pargi (nc)
+		call pargi (nl)
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+	    call error (0, Memc[errstr])
+	}
+	if ((c1 == 1) && (c2 == nc) && (l1 == 1) && (l2 == nl)) {
+	    call error (0, "Bias section not specified or given as full image")
+	}
+
+	# If no processing is desired then print overscan strip and return.
+	if (clgetb ("noproc")) {
+	    call eprintf ("  [TO BE DONE] Overscan section is [%d:%d,%d:%d].\n")
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+		call sfree (sp)
+		return
+	}
+
+	# Determine the overscan section parameters. The readout axis
+	# determines the type of overscan.  The step sizes are ignored.
+	# The limits in the long dimension are replaced by the trim limits.
+
+	type = clgwrd ("function", Memc[func], SZ_LINE, OVERSCAN_TYPES)
+	if (type < OVERSCAN_FIT) {
+	    overscan = NULL
+	    if (READAXIS(ccd) == 2)
+		call error (1,
+		    "Overscan function type not allowed with readaxis of 2")
+	} else {
+	    if (READAXIS(ccd) == 1) {
+		first = c1
+		last = c2
+		navg = last - first + 1
+		npts = nl
+		call salloc (buf, npts, TY_REAL)
+		do i = 1, npts
+		    Memr[buf+i-1] = asumr (Memr[imgs2r (IN_IM(ccd), first, last,
+			i, i)], navg)
+		if (navg > 1)
+		    call adivkr (Memr[buf], real (navg), Memr[buf], npts)
+
+		# Trim the overscan vector and set the pixel coordinate.
+		npts = CCD_L2(ccd) - CCD_L1(ccd) + 1
+		call malloc (overscan, npts, TY_REAL)
+		call salloc (x, npts, TY_REAL)
+		call trim_overscan (Memr[buf], npts, IN_L1(ccd), Memr[x],
+		    Memr[overscan])
+
+		call fit_overscan (Memc[str], c1, c2, l1, l2, Memr[x],
+		    Memr[overscan], npts)
+
+	    } else {
+		first = l1
+		last = l2
+		navg = last - first + 1
+		npts = nc
+		call salloc (buf, npts, TY_REAL)
+		call aclrr (Memr[buf], npts)
+		do i = first, last
+		    call aaddr (Memr[imgl2r(IN_IM(ccd),i)], Memr[buf],
+			Memr[buf], npts)
+		if (navg > 1)
+		    call adivkr (Memr[buf], real (navg), Memr[buf], npts)
+
+		# Trim the overscan vector and set the pixel coordinate.
+		npts = CCD_C2(ccd) - CCD_C1(ccd) + 1
+		call malloc (overscan, npts, TY_REAL)
+		call salloc (x, npts, TY_REAL)
+		call trim_overscan (Memr[buf], npts, IN_C1(ccd), Memr[x],
+		    Memr[overscan])
+
+		call fit_overscan (Memc[str], c1, c2, l1, l2, Memr[x],
+		    Memr[overscan], npts)
+	    }
+	}
+	
+	# Set the CCD structure overscan parameters.
+	CORS(ccd, OVERSCAN) = O
+	COR(ccd) = YES
+	OVERSCAN_TYPE(ccd) = type
+	OVERSCAN_VEC(ccd) = overscan
+
+	# Log the operation.
+	if (type < OVERSCAN_FIT) {
+	    call sprintf (Memc[str], SZ_LINE,
+		"Overscan section is [%d:%d,%d:%d] with function=%s")
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+		call pargstr (Memc[func])
+	} else {
+	    call sprintf (Memc[str], SZ_LINE,
+		"Overscan section is [%d:%d,%d:%d] with mean=%g")
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+		call pargr (asumr (Memr[overscan], npts) / npts)
+	}
+	call timelog (Memc[str], SZ_LINE)
+	call ccdlog (IN_IM(ccd), Memc[str])
+	call hdmpstr (OUT_IM(ccd), "overscan", Memc[str])
+
+	call sfree (sp)
+end
+
+
+# FIT_OVERSCAN -- Fit a function to smooth the overscan vector.
+#   The fitting uses the ICFIT procedures which may be interactive.
+#   Changes to these parameters are "learned".  The user is queried with a four
+#   valued logical query (XT_ANSWER routine) which may be turned off when
+#   multiple images are processed.
+
+procedure fit_overscan (image, c1, c2, l1, l2, x, overscan, npts)
+
+char	image[ARB]		# Image name for query and title
+int	c1, c2, l1, l2		# Overscan strip
+real	x[npts]			# Pixel coordinates of overscan
+real	overscan[npts]		# Input overscan and output fitted overscan
+int	npts			# Number of data points
+
+int	interactive, fd
+pointer	sp, str, w, ic, cv, gp, gt
+
+int	clgeti(), ic_geti(), open()
+real	clgetr(), ic_getr()
+pointer	gopen(), gt_init()
+errchk	gopen, open
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (w, npts, TY_REAL)
+	call amovkr (1., Memr[w], npts)
+
+	# Open the ICFIT procedures, get the fitting parameters, and
+	# set the fitting limits.
+
+	call ic_open (ic)
+	call clgstr ("function", Memc[str], SZ_LINE)
+	call ic_pstr (ic, "function", Memc[str])
+	call ic_puti (ic, "order", clgeti ("order"))
+	call clgstr ("sample", Memc[str], SZ_LINE)
+	call ic_pstr (ic, "sample", Memc[str])
+	call ic_puti (ic, "naverage", clgeti ("naverage"))
+	call ic_puti (ic, "niterate", clgeti ("niterate"))
+	call ic_putr (ic, "low", clgetr ("low_reject"))
+	call ic_putr (ic, "high", clgetr ("high_reject"))
+	call ic_putr (ic, "grow", clgetr ("grow"))
+	call ic_putr (ic, "xmin", min (x[1], x[npts]))
+	call ic_putr (ic, "xmax", max (x[1], x[npts]))
+	call ic_pstr (ic, "xlabel", "Pixel")
+	call ic_pstr (ic, "ylabel", "Overscan")
+
+	# If the fitting is done interactively set the GTOOLS and GIO
+	# pointers.  Also "learn" the fitting parameters since they may
+	# be changed when fitting interactively.
+
+	call sprintf (Memc[str], SZ_LINE,
+	    "Fit overscan vector for %s interactively")
+	    call pargstr (image)
+	call set_interactive (Memc[str], interactive)
+	if ((interactive == YES) || (interactive == ALWAYSYES)) {
+	    gt = gt_init ()
+	    call sprintf (Memc[str], SZ_LINE,
+	        "Overscan vector for %s from section [%d:%d,%d:%d]\n")
+	        call pargstr (image)
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+	    call gt_sets (gt, GTTITLE, Memc[str])
+	    call gt_sets (gt, GTTYPE, "line")
+	    call gt_setr (gt, GTXMIN, x[1])
+	    call gt_setr (gt, GTXMAX, x[npts])
+	    call clgstr ("graphics", Memc[str], SZ_FNAME)
+	    gp = gopen (Memc[str], NEW_FILE, STDGRAPH)
+
+	    call icg_fit (ic, gp, "cursor", gt, cv, x, overscan, Memr[w], npts)
+
+	    call ic_gstr (ic, "function", Memc[str], SZ_LINE)
+	    call clpstr ("function", Memc[str])
+	    call clputi ("order", ic_geti (ic, "order"))
+	    call ic_gstr (ic, "sample", Memc[str], SZ_LINE)
+	    call clpstr ("sample", Memc[str])
+	    call clputi ("naverage", ic_geti (ic, "naverage"))
+	    call clputi ("niterate", ic_geti (ic, "niterate"))
+	    call clputr ("low_reject", ic_getr (ic, "low"))
+	    call clputr ("high_reject", ic_getr (ic, "high"))
+	    call clputr ("grow", ic_getr (ic, "grow"))
+
+	    call gclose (gp)
+	    call gt_free (gt)
+	} else
+	    call ic_fit (ic, cv, x, overscan, Memr[w], npts, YES, YES, YES, YES)
+
+	# Make a log of the fit in the plot file if given.
+	call clgstr ("plotfile", Memc[str], SZ_LINE)
+	call xt_stripwhite (Memc[str])
+	if (Memc[str] != EOS) {
+	    fd = open (Memc[str], APPEND, BINARY_FILE)
+	    gp = gopen ("stdvdm", NEW_FILE, fd)
+	    gt = gt_init ()
+	    call sprintf (Memc[str], SZ_LINE,
+	        "Overscan vector for %s from section [%d:%d,%d:%d]\n")
+	        call pargstr (image)
+		call pargi (c1)
+		call pargi (c2)
+		call pargi (l1)
+		call pargi (l2)
+	    call gt_sets (gt, GTTITLE, Memc[str])
+	    call gt_sets (gt, GTTYPE, "line")
+	    call gt_setr (gt, GTXMIN, 1.)
+	    call gt_setr (gt, GTXMAX, real (npts))
+	    call icg_graphr (ic, gp, gt, cv, x, overscan, Memr[w], npts)
+	    call gclose (gp)
+	    call close (fd)
+	    call gt_free (gt)
+	}
+
+	# Replace the raw overscan vector with the smooth fit.
+	call cvvector (cv, x, overscan, npts)
+
+	# Finish up.
+	call ic_closer (ic)
+	call cvfree (cv)
+	call sfree (sp)
+end
+
+
+# TRIM_OVERSCAN -- Trim the overscan vector.
+
+procedure trim_overscan (data, npts, start, x, overscan)
+
+real	data[ARB]		# Full overscan vector
+int	npts			# Length of trimmed vector
+int	start			# Trim start
+real	x[npts]			# Trimmed pixel coordinates (returned)
+real	overscan[npts]		# Trimmed overscan vector (returned)
+
+int	i, j
+
+begin
+	do i = 1, npts {
+	    j = start + i - 1
+	    x[i] = j
+	    overscan[i] = data[j]
+	}
+end