Repatch (from linux) of OSX IRAF

1 files changed, 344 insertions, 0 deletions

diff --git a/noao/imred/quadred/src/ccdproc/setoverscan.x b/noao/imred/quadred/src/ccdproc/setoverscan.x
new file mode 100644
index 00000000..2fef378a
--- /dev/null
+++ b/noao/imred/quadred/src/ccdproc/setoverscan.x
@@ -0,0 +1,344 @@
+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.	Average the overscan columns or lines.
+#   4.	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, j, nsec, navg, npts, first, last
+int	nc, nl, c1, c2, l1, l2
+pointer	sp, str, errstr, buf, overscan, x, y, z
+
+real	asumr()
+bool	clgetb(), ccdflag()
+pointer	imgl2r()
+errchk	imgl2r, 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 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)
+	navg = c2 - c1 + 1
+	npts = CCD_L2(ccd) - CCD_L1(ccd) + 1
+
+	nsec = max (1, IN_NSEC(ccd))
+	do i = 1, nsec {
+	    if (BIAS_SEC(ccd) != NULL) {
+		c1 = BIAS_SC1(ccd,i)
+		c2 = BIAS_SC2(ccd,i)
+		l1 = BIAS_SL1(ccd,i)
+		l2 = BIAS_SL2(ccd,i)
+	    }
+	    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.
+
+	if (READAXIS(ccd) == 1) {
+	    call salloc (buf, nsec*nl, TY_REAL)
+	    z = buf
+	    do i = 1, nl {
+		y = imgl2r (IN_IM(ccd), i)
+		do j = 1, nsec {
+		    if (BIAS_SEC(ccd) != NULL) {
+			l1 = BIAS_SL1(ccd,j)
+			l2 = BIAS_SL2(ccd,j)
+			if (i < l1 || i > l2)
+			    next
+			c1 = BIAS_SC1(ccd,j)
+			c2 = BIAS_SC2(ccd,j)
+			navg = c2 - c1 + 1
+			z = buf + (j - 1) * nl
+		    }
+		    Memr[z+i-1] = asumr (Memr[y+c1-1], navg)
+		}
+	    }
+
+	    # Trim the overscan vector and set the pixel coordinate.
+	    call salloc (x, nl, TY_REAL)
+	    call malloc (overscan, nsec*nl, TY_REAL)
+	    y = overscan
+	    do i = 1, nsec {
+		if (BIAS_SEC(ccd) != NULL) {
+		    c1 = BIAS_SC1(ccd,i)
+		    c2 = BIAS_SC2(ccd,i)
+		    l1 = BIAS_SL1(ccd,i)
+		    l2 = BIAS_SL2(ccd,i)
+		    navg = c2 - c1 + 1
+		    npts = l2 - l1 + 1
+		    y = overscan + (i - 1) * nl
+		    z = buf + (i - 1) * nl
+		}
+		if (navg > 1)
+		    call adivkr (Memr[z+l1-1], real (navg), Memr[z+l1-1],
+			npts)
+		call trim_overscan (Memr[z], npts, l1, Memr[x], Memr[y])
+		call fit_overscan (Memc[str], c1, c2, l1, l2, Memr[x],
+		    Memr[y], 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_VEC(ccd) = overscan
+
+	# Log the operation.
+	call strcpy ("overscan", Memc[errstr], SZ_LINE)
+	y = overscan
+	do i = 1, nsec {
+	    if (BIAS_SEC(ccd) != NULL) {
+		c1 = BIAS_SC1(ccd,i)
+		c2 = BIAS_SC2(ccd,i)
+		l1 = BIAS_SL1(ccd,i)
+		l2 = BIAS_SL2(ccd,i)
+		y = overscan + (i - 1) * nl
+		npts = c2 - c1 + 1
+		if (i > 1) {
+		    call sprintf (Memc[errstr], SZ_LINE, "ovrscn%d")
+			call pargi (i)
+		}
+	    }
+	    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[y], npts) / npts)
+	    call timelog (Memc[str], SZ_LINE)
+	    call ccdlog (IN_IM(ccd), Memc[str])
+	    call hdmpstr (OUT_IM(ccd), Memc[errstr], 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[%d:%d,%d:%d] interactively")
+	    call pargstr (image)
+	    call pargi (c1)
+	    call pargi (c2)
+	    call pargi (l1)
+	    call pargi (l2)
+	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