1 files changed, 229 insertions, 0 deletions
diff --git a/noao/imred/ccdred/src/t_mkillumft.x b/noao/imred/ccdred/src/t_mkillumft.x
new file mode 100644
index 00000000..ecb66a8e
--- /dev/null
+++ b/noao/imred/ccdred/src/t_mkillumft.x
@@ -0,0 +1,229 @@
+include <imhdr.h>
+include	"ccdred.h"
+
+
+# T_MKILLUMFLAT -- Make illumination corrected flat field images.
+#
+# The input flat field images are processed and smoothed to obtain
+# illumination pattern.  The illumination pattern is then divided out
+# of the input image to make the output illumination corrected flat field
+# image.
+
+procedure t_mkillumflat()
+
+int	listin			# List of input CCD images
+int	listout			# List of output CCD images
+int	ccdtype			# CCD image type
+int	interactive		# Fit overscan interactively?
+
+bool	clgetb(), streq()
+int	imtopenp(), imtgetim()
+pointer	sp, input, output, tmp, str, in, out, ccd
+errchk	set_input, set_output, ccddelete
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (tmp, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the lists and instrument translation file.  Open the translation
+	# file.  Initialize the interactive flag and the calibration images.
+
+	listin = imtopenp ("input")
+	listout = imtopenp ("mkillumflat.output")
+	call clgstr ("instrument", Memc[input], SZ_FNAME)
+	call hdmopen (Memc[input])
+	call set_interactive ("", interactive)
+	call cal_open (NULL)
+	call ccd_open (0)
+
+	# Process each image.
+	while (imtgetim (listin, Memc[input], SZ_FNAME) != EOF) {
+	    if (clgetb ("noproc")) {
+		call printf ("%s: mkillumflat\n")
+		    call pargstr (Memc[input])
+	    }
+	    
+	    # Set input and output images.  Use temporary image if needed.
+	    call set_input (Memc[input], in, ccdtype)
+	    if (in == NULL)
+		next
+
+	    if (imtgetim (listout, Memc[output], SZ_FNAME) == EOF)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+	    if (Memc[output] == EOS)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+	    if (streq (Memc[input], Memc[output]))
+	        call mktemp ("tmp", Memc[tmp], SZ_FNAME)
+	    else
+		call strcpy (Memc[output], Memc[tmp], SZ_FNAME)
+	    call set_output (in, out, Memc[tmp])
+
+	    # Process image as a flat field image.
+	    call set_proc (in, out, ccd)
+	    call set_sections (ccd)
+	    call set_trim (ccd)
+	    call set_fixpix (ccd)
+	    call set_overscan (ccd)
+	    call set_zero (ccd)
+	    call set_dark (ccd)
+
+	    # Do the processing.
+	    if (CORS(ccd) == YES) {
+	        call doproc (ccd)
+	        call set_header (ccd)
+
+	        # Finish up
+	        call imunmap (in)
+	        call imunmap (out)
+	        if (streq (Memc[input], Memc[output])) {
+		    call ccddelete (Memc[input])
+	            call imrename (Memc[tmp], Memc[input])
+		} else
+		    call strcpy (Memc[output], Memc[input], SZ_FNAME)
+	    } else {
+		# Delete the temporary output image.  Make a copy if needed.
+	        call imunmap (in)
+	        call imunmap (out)
+		call imdelete (Memc[tmp])
+	    }
+	    call free_proc (ccd)
+
+	    # Do special processing.
+	    call mkillumflat (Memc[input], Memc[output])
+	    if (!streq (Memc[input], Memc[output]))
+		call ccdcopy (Memc[input], Memc[output])
+	}
+
+	# Finish up.
+	call hdmclose ()
+	call imtclose (listin)
+	call imtclose (listout)
+	call cal_close ()
+	call ccd_close ()
+	call sfree (sp)
+end
+
+
+# MKILLUMFLAT -- Take the processed input image and make the illumination
+# corrected flat field output image.  The illumination pattern is created
+# as a temporary image and then the applied to the input flat field
+# image to make the final output flat field image.  If the input and
+# output names are the same the operation is done in place.
+
+procedure mkillumflat (input, output)
+
+char	input[SZ_FNAME]		# Input image
+char	output[SZ_FNAME]	# Output image
+
+int	i, nc, nl
+real	scale
+long	time
+pointer	sp, str, illum, tmp, in, im, out, out1, data
+
+bool	clgetb(), ccdflag(), streq()
+int	hdmgeti()
+real	hdmgetr(), clgetr(), divzero()
+pointer	immap(), imgl2r(), impl2r()
+errchk	immap, ccddelete
+extern	divzero()
+
+real	rdivzero	# Result for divion by zero
+int	ndivzero	# Number of zero divisions
+common	/cdivzero/ rdivzero, ndivzero
+
+begin
+	# Check if this operation has been done.
+	in = immap (input, READ_ONLY, 0)
+	if (ccdflag (in, "illumflt")) {
+	    call imunmap (in)
+	    return
+	}
+
+	# Print operation if not processing.
+	if (clgetb ("noproc")) {
+	    call eprintf (
+		"  [TO BE DONE] Remove illumination\n")
+		call pargstr (input)
+	    call imunmap (in)
+	    return
+	}
+
+	# Get and set task parameters for division by zero.
+	rdivzero = clgetr ("divbyzero")
+	ndivzero = 0
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (illum, SZ_FNAME, TY_CHAR)
+	call salloc (tmp, SZ_FNAME, TY_CHAR)
+
+	# Make the illumination image.
+	call imunmap (in)
+	call strcpy (input, Memc[tmp], SZ_FNAME)
+	call mktemp ("tmp", Memc[illum], SZ_FNAME)
+	call mkillumination (Memc[tmp], Memc[illum], NO, NO)
+
+	in = immap (input, READ_ONLY, 0)
+	im = immap (Memc[illum], READ_ONLY, 0)
+	iferr (scale = hdmgetr (im, "ccdmean"))
+	    scale = 1.
+	iferr (time = hdmgeti (im, "ccdmeant"))
+	    time = IM_MTIME(im)
+	if (time < IM_MTIME(im))
+	    scale = 1.
+
+	# Create the temporary output.
+	if (streq (input, output)) {
+	    call mktemp ("tmp", Memc[tmp], SZ_FNAME)
+	    call set_output (in, out, Memc[tmp])
+	    out1 = in
+	} else {
+	    call set_output (in, out, output)
+	    out1 = out
+	}
+
+	# Divide the illumination and flat field images with scaling.
+	nc = IM_LEN(out,1)
+	nl = IM_LEN(out,2)
+	do i = 1, nl {
+	    data = impl2r (out, i)
+	    call advzr (Memr[imgl2r(in,i)], Memr[imgl2r(im,i)],
+		Memr[data], nc, divzero)
+	    if (scale != 1.)
+	        call amulkr (Memr[data], scale, Memr[data], nc)
+	}
+
+	# Log the operation.
+	if (ndivzero > 0) {
+	    call sprintf (Memc[str], SZ_LINE,
+		"Warning: %d divisions by zero replaced by %g")
+		call pargi (ndivzero)
+		call pargr (rdivzero)
+	    call ccdlog (out1, Memc[str])
+	}
+	call sprintf (Memc[str], SZ_LINE, "Removed illumination from flat")
+	call sprintf (Memc[str], SZ_LINE,
+	    "Illumination flat created from %s")
+	    call pargstr (input)
+	call timelog (Memc[str], SZ_LINE)
+	call ccdlog (out1, Memc[str])
+	call hdmpstr (out, "illumflt", Memc[str])
+	call hdmpstr (out, "imagetyp", "flat")
+
+	# Finish up
+	call imunmap (in)
+	call imunmap (im)
+	call imunmap (out)
+	call imdelete (Memc[illum])
+
+	# The input name is changed to the output name for further processing.
+	if (streq (input, output)) {
+	    call ccddelete (input)
+	    call imrename (Memc[tmp], input)
+	} else
+	    call strcpy (output, input, SZ_FNAME)
+	call sfree (sp)
+end