Repatch (from linux) of OSX IRAF

8 files changed, 1922 insertions, 0 deletions

diff --git a/noao/imred/ccdred/src/cosmic/cosmicrays.hlp b/noao/imred/ccdred/src/cosmic/cosmicrays.hlp
new file mode 100644
index 00000000..bfb56e9c
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/cosmicrays.hlp
@@ -0,0 +1,338 @@
+.help cosmicrays Dec87 noao.imred.ccdred
+.ih
+NAME
+cosmicrays -- Detect and replace cosmic rays
+.ih
+USAGE
+cosmicrays input output
+.ih
+PARAMETERS
+.ls input
+List of input images in which to detect cosmic rays.
+.le
+.ls output
+List of output images in which the detected cosmic rays will be replaced
+by an average of neighboring pixels.  If the output image name differs
+from the input image name then a copy of the input image is made with
+the detected cosmic rays replaced.  If no output images are specified
+then the input images are modified in place.  In place modification of
+an input image also occurs when the output image name is the same as
+the input image name.
+.le
+.ls badpix = ""
+List of bad pixel files to be created, one for each input image.  If no
+file names are given then no bad pixel file is created.  The bad pixel 
+file is a simple list of pixel coordinates for each replaced cosmic ray.
+This file may be used in conjunction with \fBbadpixelimage\fR to create
+a mask image.
+.le
+
+.ls ccdtype = ""
+If specified only the input images of the desired CCD image type will be
+selected.
+.le
+.ls threshold = 25.
+Detection threshold above the mean of the surrounding pixels for cosmic
+rays.  The threshold will depend on the noise characteristics of the
+image and how weak the cosmic rays may be for detection.  A typical value
+is 5 or more times the sigma of the background.
+.le
+.ls fluxratio = 2.
+The ratio (as a percent) of the mean neighboring pixel flux to the candidate
+cosmic ray pixel for rejection.  The value depends on the seeing and the
+characteristics of the cosmic rays.  Typical values are in the range
+2 to 10 percent.  This value may be reset interactively from a plot
+or defined by identifying selected objects as stars or cosmic rays.
+.le
+.ls npasses = 5
+Number of cosmic ray detection passes.  Since only the locally strongest
+pixel is considered a cosmic ray, multiple detection passes are needed to
+detect and replace multiple pixel cosmic ray events.
+.le
+.ls window = 5
+Size of cosmic ray detection window.  A square window of either 5 by 5 or
+7 by 7 is used to detect cosmic rays.  The smaller window allows detection
+in the presence of greater background gradients but is less sensitive at
+discriminating multiple event cosmic rays from stars.  It is also marginally
+faster.
+.le
+.ls interactive = yes
+Examine parameters interactively?  A plot of the mean flux within the
+detection window (x100) vs the flux ratio (x100) is plotted and the user may
+set the flux ratio threshold, delete and undelete specific events, and
+examine specific events.  This is useful for new data in which one is
+uncertain of an appropriate flux ratio threshold.  Once determined the
+task need not be used interactively.
+.le
+.ls train = no
+Define the flux ratio threshold by using a set of objects identified
+as stars (or other astronomical objects) or cosmic rays?
+.le
+.ls objects = ""
+Cursor list of coordinates of training objects.  If null (the null string "")
+then the image display cursor will be read.  The user is responsible for first
+displaying the image.  Otherwise a file containing cursor coordinates
+may be given.  The format of the cursor file is "x y wcs key" where
+x and y are the pixel coordinates, wcs is an arbitrary number such as 1,
+and key may be 's' for star or 'c' for cosmic ray.
+.le
+.ls savefile = ""
+File to save (by appending) the training object coordinates.  This is of
+use when the objects are identified using the image display cursor.  The
+saved file can then be input as the object cursor list for repeating the
+execution.
+.le
+.ls answer
+This parameter is used for interactive queries when processing a list of
+images.  The responses may be "no", "yes", "NO", or "YES".  The upper case
+responses permanently enable or disable the interactive review while
+the lower case reponses allow selective examination of certain input
+images.  \fIThis parameter should not be specified on the command line.
+If it is then the value will be ignored and the task will act as if
+the answer "yes" is given for each image; i.e. it will enter the interactive
+phase without prompting.\fR
+.le
+.ih
+OTHER PARAMETERS
+There are other parameters which may be defined by the package, as is the
+case with \fBccdred\fR, or as part of the task, as is the case with
+standalone version in the \fBgeneric\fR package.
+
+.ls verbose
+If yes then a time stamped log of the operation is printed on the standard
+output.
+.le
+.ls logfile
+If a log file is specified then a time stamped log of the operation is
+recorded.
+.le
+.ls plotfile
+If a plot file is specified then the graph of the flux ratio (x100) vs
+the mean flux (x100) is recorded as metacode.  This may be spooled or examined
+later.
+.le
+.ls graphics = "stdgraph"
+Interactive graphic output device for interactive examination of the
+detection parameters.
+.le
+.ls cursor = ""
+Interactive graphics cursor input.  If null the graphics display cursor
+is used, otherwise a file containing cursor input may be specified.
+.le
+.ls instrument
+The \fBccdred\fR instrument file is used for mapping header keywords and
+CCD image types.
+.le
+.ih
+IMAGE CURSOR COMMANDS
+
+.nf
+?	Help
+c	Identify the object as a cosmic ray
+s	Identify the object as a star
+g	Switch to the graphics plot
+q	Quit and continue with the cleaning
+.fi
+
+GRAPHICS CURSOR COMMANDS
+
+.nf
+?	Help
+a	Toggle between showing all candidates and only the training points
+d	Mark candidate for replacement (applys to '+' points)
+q	Quit and return to image cursor or replace the selected pixels
+r	Redraw the graph
+s	Make a surface plot for the candidate nearest the cursor
+t	Set the flux ratio threshold at the y cursor position
+u	Mark candidate to not be replaced (applys to 'x' points)
+w	Adjust the graph window (see \fBgtools\fR)
+<space>	Print the pixel coordinates
+.fi
+
+There are no colon commands except those for the windowing options (type
+:\help or see \fBgtools\fR).
+.ih
+DESCRIPTION
+Cosmic ray events in each input image are detected and replaced by the
+average of the four neighbors.  The replacement may be performed
+directly on the input image if no output image is specified or if the
+output image name is the same as the input image name.  If a new image
+is created it is a copy of the input image except for the replaced
+pixels.  The processing keyword CRCOR is added to the output image
+header.  Optional output includes a log file to which a processing log
+is appended, a verbose log output to the standard output (the same as
+that in the log file), a plot file showing the parameters of the
+detected cosmic ray candidates and the flux ratio threshold used, a
+bad pixel file containing the coordinates of the replaced pixels, and
+a file of training objects marked with the image display cursor.  The
+bad pixel file may be used for plotting purposes or to create a mask
+image for display and analysis using the task \fBbadpiximage\fR.  This
+bad pixel file will be replaced by the IRAF bad pixel facility when it
+becomes available.  If one wants more than a simple mask image then by
+creating a different output image a difference image between the
+original and the modified image may be made using \fBimarith\fR.
+
+This task may be applied to an image previously processed to detect
+additional cosmic rays.  A warning will be given (because of the
+CRCOR header parameter) and the previous processing header keyword will
+be overwritten.
+
+The cosmic ray detection algorithm consists of the following steps.
+First a pixel must be the brightest pixel within the specified
+detection window (either 5x5 or 7x7).  The mean flux in the surrounding
+pixels with the second brightest pixel excluded (which may also be a
+cosmic ray event) is computed and the candidate pixel must exceed this
+mean by the amount specified by the parameter \fIthreshold\fR.  A plane
+is fit to the border pixels of the window and the fitted background is
+subtracted.  The mean flux (now background subtracted) and the ratio of
+this mean to the cosmic ray candidate (the brightest pixel) are
+computed.  The mean flux (x100) and the ratio (x100) are recorded for
+interactive examination if desired.
+
+Once the list of cosmic ray candidates has been created and a threshold for
+the flux ratio established (by the parameter \fIfluxratio\fR, by the
+"training" method, or by using the graphics cursor in the interactive plot)
+the pixels with ratios below the threshold are replaced in the output by
+the average of the four neighboring pixels (with the second strongest pixel
+in the detection window excluded if it is one of these pixels).  Additonal
+pixels may then be detected and replaced in further passes as specified by
+the parameter \fInpasses\fR.  Note that only pixels in the vicinity of
+replaced pixels need be considered in further passes.
+
+The division between the peaks of real objects and cosmic rays is made
+based on the flux ratio between the mean flux (excluding the center
+pixel and the second strongest pixel) and the candidate pixel.  This
+threshold depends on the point spread function and the distribution of
+multiple cosmic ray events and any additional neighboring light caused
+by the events.  This threshold is not strongly coupled to small changes
+in the data so that once it is set for a new type of image data it may
+be used for similar images.  To set it initially one may examine the
+scatter plot of the flux ratio as a function of the mean flux.  This
+may be done interactively or from the optional plot file produced.
+
+After the initial list of cosmic ray candidates has been created and before
+the final replacing cosmic rays there are two optional steps to allow
+examining the candidates and setting the flux ratio threshold dividing
+cosmic rays from real objects.  The first optional step is define the flux
+ratio boundary by reference to user specified classifications; that is
+"training".  To do this step the \fItrain\fR parameter must be set to yes.
+The user classified objects are specified by a cursor input list.  This
+list can be an actual file or the image display cursor as defined by the
+\fIobjects\fR parameter.  The \fIsavefile\fR parameter is also used during
+the training to record the objects specified.  The parameter specifies a
+file to append the objects selected.  This is useful when the objects are
+defined by interactive image cursor and does not make much sense when using
+an input list.
+
+If the \fIobjects\fR parameter is specified as a null string then
+the image display cursor will be repeatedly read until a 'q' is
+entered.  The user first displays the image and then when the task
+reads the display cursor the cursor shape will change.  The user
+points at objects and types 's' for a star (or other astronomical
+object) and 'c' for a cosmic ray.  Note that this input is used
+to search for the matching object in the cosmic ray candidate list
+and so it is possible the selected object is not in the list though
+it is unlikely.  The selection will be quietly ignored in that case.
+To exit the interactive selection of training objects type 'q'.
+
+If 'g' is typed a graph of all the candidates is drawn showing
+"flux" vs. "flux ratio" (see below for more).  Training objects will
+be shown with a box and the currently set flux ratio threshold will
+also be shown.  Exiting the plot will return to entering more training
+objects.  The plot will remain and additional objects will immediately
+be shown with a new box.  Thus, if one wants to see the training
+objects identified in the plot as one selects them from the image
+display first type a 'g' to draw the initial plot.  Also by switching
+to the plot with 'g' allows you to draw surface plots (with 's') or
+get the pixel coordinates of a candidate (the space key) to be
+found in the display using the coordinate readout of the display.
+Note that the display interaction is simpler than might be desired
+because this task does not directly connect to the display.
+
+The most likely use for training is with the interactive image display.
+However one may prepare an input list by other means, one example
+is with \fBrimcursor\fR, and then specify the file name.  The savefile
+may also be used a cursor input to repeat the cosmic ray operation
+(but be careful not to have the cursor input and save file be the
+same file!).
+
+The flux ratio threshold is determined from the training objects by
+finding the point with the minimum number of misclassifications
+(stars as cosmic rays or cosmic rays as stars).  The threshold is
+set at the lowest value so that it will always go through one of
+the cosmic ray objects.  There should be at least one of each type
+of object defined for this to work.  The following option of
+examining the cosmic ray candidates and parameters may still be
+used to modify the derived flux ratio threshold.  One last point
+about the training objects is that even if some of the points
+lie on the wrong side of the threshold they will remain classified
+as cosmic ray or non-cosmic ray.  In other words, any object
+classified by the user will remain in that classification regardless
+of the final flux ratio threshold.
+
+After the training step the user will be queried to examine the candidates
+in the flux vs flux ratio plane if the \fIinteractive\fR flag is set.
+Responses may be made for specific images or for all images by using
+lower or upper case answers respectively.  When the parameters are
+examined interactively the user may change the flux ratio threshold
+('t' key).  Changes made are stored in the parameter file and, thus,
+learned for further images.  Pixels to be deleted are marked by crosses
+and pixels which are peaks of objects are marked by pluses.  The user
+may explicitly delete or undelete any point if desired but this is only
+for special cases near the threshold.  In the future keys for
+interactive display of the specific detections will be added.
+Currently a surface plot of any candidate may be displayed graphically
+in four 90 degree rotated views using the 's' key.  Note that the
+initial graph does not show all the points some of which are clearly
+cosmic rays because they have negative mean flux or flux ratio.  To
+view all data one must rewindow the graph with the 'w' key or ":/"
+commands (see \fBgtools\fR).
+.ih
+EXAMPLES
+1. To replace cosmic rays in a set of images ccd* without training:
+
+.nf
+    cl> cosmicrays ccd* new//ccd*
+    ccd001: Examine parameters interactively? (yes):
+    [A scatter plot graph is made.  One can adjust the threshold.]
+    [Looking at a few points using the 's' key can be instructive.]
+    [When done type 'q'.]
+    ccd002: Examine parameters interactively? (yes): NO
+    [No further interactive examination is done.]
+.fi
+
+After cleaning one typically displays the images and  possibly blinks them.
+A difference image or mask image may also be created.
+
+2. To use the interactive training method for setting the flux ratio threshold:
+
+.nf
+    # First display the image.
+    cl> display ccd001 1
+    z1 = 123.45 z2= 543.21
+    cl> cosmicrays ccd001 ccd001cr train+
+    [After the cosmic ray candidates are found the image display
+    [cursor will be activated.  Mark a cosmic ray with 'c' and
+    [a star with 's'.  Type 'g' to get a plot showing the two
+    [points with boxes.  Type 'q' to go back to the image display.
+    [As each new object is marked a box will appear in the plot and
+    [the threshold may change.  To find the location of an object
+    [seen in the plot use 'g' to go to the graph, space key to find
+    [the pixel coordinates, 'q' to go back to the image display,
+    [and the image display coordinate box to find the object.
+    [When done with the training type 'q'.
+    ccd001: Examine parameters interactively? (yes): no
+.fi
+
+3.  To create a mask image a bad pixel file must be specified.  In the
+following we replace the cosmic rays in place and create a bad pixel
+file and mask image:
+
+.nf
+    cl> cosmicrays ccd001 ccd001 badpix=ccd001.bp
+    cl> badpiximage ccd001.bp ccd001 ccd001bp
+.fi
+.ih
+SEE ALSO
+badpixelimage gtools imedit rimcursor
+.endhelp
diff --git a/noao/imred/ccdred/src/cosmic/crexamine.x b/noao/imred/ccdred/src/cosmic/crexamine.x
new file mode 100644
index 00000000..d84961bc
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/crexamine.x
@@ -0,0 +1,486 @@
+include	<error.h>
+include	<syserr.h>
+include <imhdr.h>
+include	<gset.h>
+include	<mach.h>
+include	<pkg/gtools.h>
+include	"crlist.h"
+
+# CR_EXAMINE  -- Examine cosmic ray candidates interactively.
+# CR_GRAPH    -- Make a graph
+# CR_NEAREST  -- Find the nearest cosmic ray to the cursor.
+# CR_DELETE   -- Set replace flag for cosmic ray candidate nearest cursor.
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+# CR_UPDATE   -- Change replacement flags, thresholds, and graphs.
+# CR_PLOT     -- Make log plot
+
+define	HELP	"noao$lib/scr/cosmicrays.key"
+define	PROMPT	"cosmic ray options"
+
+# CR_EXAMINE -- Examine cosmic ray candidates interactively.
+
+procedure cr_examine (cr, gp, gt, im, fluxratio, first)
+
+pointer	cr			# Cosmic ray list
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+pointer	im			# Image pointer
+real	fluxratio		# Flux ratio threshold
+int	first			# Initial key
+
+char	cmd[SZ_LINE]
+int	i, newgraph, wcs, key, nc, nl, c1, c2, l1, l2, show
+real	wx, wy
+pointer	data
+
+int	clgcur()
+pointer	imgs2r()
+
+begin
+	# Set up the graphics.
+	call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+
+	# Set image limits
+	nc = IM_LEN(im, 1)
+	nl = IM_LEN(im, 2)
+
+	# Enter cursor loop.
+	key = first
+	repeat {
+	    switch (key) {
+	    case '?': # Print help text.
+		call gpagefile (gp, HELP, PROMPT)
+	    case ':': # Colon commands.
+		switch (cmd[1]) {
+		case '/':
+		    call gt_colon (cmd, gp, gt, newgraph)
+		default:
+		    call printf ("\007")
+		}
+	    case 'a': # Toggle show all
+		if (show == 0)
+		    show = 1
+		else
+		    show = 0
+		newgraph = YES
+	    case 'd': # Delete candidate
+		call cr_delete (gp, wx, wy, cr, i, show)
+	    case 'q': # Quit
+		break
+	    case 'r': # Redraw the graph.
+		newgraph = YES
+	    case 's': # Make surface plots
+		call cr_nearest (gp, wx, wy, cr, i, show)
+		c1 = max (1, int (Memr[CR_COL(cr)+i-1]) - 5)
+		c2 = min (nc, int (Memr[CR_COL(cr)+i-1]) + 5)
+		l1 = max (1, int (Memr[CR_LINE(cr)+i-1]) - 5)
+		l2 = min (nl, int (Memr[CR_LINE(cr)+i-1]) + 5)
+		data = imgs2r (im, c1, c2, l1, l2)
+		call gclear (gp)
+		call gsview (gp, 0.03, 0.48, 0.53, 0.98)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -33., 25.)
+		call gsview (gp, 0.53, 0.98, 0.53, 0.98)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -123., 25.)
+		call gsview (gp, 0.03, 0.48, 0.03, 0.48)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 57., 25.)
+		call gsview (gp, 0.53, 0.98, 0.03, 0.48)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 147., 25.)
+		call fprintf (STDERR, "[Type any key to continue]")
+		i = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)
+		newgraph = YES
+	    case 't': # Set threshold
+		call cr_update (gp, wy, cr, fluxratio, show)
+		call clputr ("fluxratio", fluxratio)
+	    case 'u': # Undelete candidate
+		call cr_undelete (gp, wx, wy, cr, i, show)
+	    case 'w':# Window the graph.
+		call gt_window (gt, gp, "cursor", newgraph)
+	    case ' ': # Print info
+		call cr_nearest (gp, wx, wy, cr, i, show)
+		call printf ("%d %d\n")
+		    call pargr (Memr[CR_COL(cr)+i-1])
+		    call pargr (Memr[CR_LINE(cr)+i-1])
+	    case 'z': # NOP
+		newgraph = NO
+	    default: # Ring bell for unrecognized commands.
+		call printf ("\007")
+	    }
+
+	    # Update the graph if needed.
+	    if (newgraph == YES) {
+		call cr_graph (gp, gt, cr, fluxratio, show)
+	        newgraph = NO
+	    }
+	} until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+end
+
+
+# CR_GRAPH -- Make a graph
+
+procedure cr_graph (gp, gt, cr, fluxratio, show)
+
+pointer	gp		# GIO pointer
+pointer	gt		# GTOOLS pointers
+pointer	cr		# Cosmic ray list
+real	fluxratio	# Flux ratio threshold
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x1, x2, y1, y2
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	call gclear (gp)
+	call gt_ascale (gp, gt, Memr[x+1], Memr[y+1], ncr)
+	call gt_swind (gp, gt)
+	call gt_labax (gp, gt)
+
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_PLUS, 2., 2.)
+	    else
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_CROSS, 2., 2.)
+	    if (Memr[w+i] != 0.)
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_BOX, 2., 2.)
+	}
+
+	call ggwind (gp, x1, x2, y1, y2)
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+
+	call sfree (sp)
+end
+
+
+# CR_NEAREST -- Find the nearest cosmic ray to the cursor.
+
+procedure cr_nearest (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+	if (index != NULL)
+	    nearest = Memi[index+nearest]
+
+	# Move the cursor to the selected point.
+	call gscur (gp, x2, y2)
+
+	call sfree (sp)
+end
+
+
+# CR_DELETE -- Set replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_delete (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	nearest = 0
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == YES) || (Memi[flag+i] == ALWAYSYES))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+
+	# Move the cursor to the selected point and mark the deleted point.
+	if (nearest > 0) {
+	    if (index != NULL)
+		nearest = Memi[index+nearest]
+	    Memi[CR_FLAG(cr)+nearest-1] = ALWAYSYES
+	    Memi[CR_WT(cr)+nearest-1] = -1
+	    call gscur (gp, x2, y2)
+	    call gseti (gp, G_PMLTYPE, 0)
+	    y2 = Memr[CR_RATIO(cr)+nearest-1]
+	    call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+	}
+
+	call sfree (sp)
+end
+
+
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_undelete (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	nearest = 0
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+
+	# Move the cursor to the selected point and mark the delete point.
+	if (nearest > 0) {
+	    if (index != NULL)
+		nearest = Memi[index+nearest]
+	    Memi[CR_FLAG(cr)+nearest-1] = ALWAYSNO
+	    Memi[CR_WT(cr)+nearest-1] = 1
+	    call gscur (gp, x2, y2)
+
+	    call gseti (gp, G_PMLTYPE, 0)
+	    y2 = Memr[CR_RATIO(cr)+nearest-1]
+	    call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+	}
+
+	call sfree (sp)
+end
+
+
+# CR_UPDATE -- Change replacement flags, thresholds, and graphs.
+
+procedure cr_update (gp, wy, cr, fluxratio, show)
+
+pointer	gp		# GIO pointer
+real	wy		# Y cursor position
+pointer	cr		# Cosmic ray list
+real	fluxratio	# Flux ratio threshold
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr, flag
+real	x1, x2, y1, y2
+pointer	x, y, f
+
+begin
+	call gseti (gp, G_PLTYPE, 0)
+	call ggwind (gp, x1, x2, y1, y2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+	fluxratio = wy
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+
+	if (show == 1)
+	    return
+
+	ncr = CR_NCR(cr)
+	x = CR_FLUX(cr) - 1
+	y = CR_RATIO(cr) - 1
+	f = CR_FLAG(cr) - 1
+
+	do i = 1, ncr {
+	    flag = Memi[f+i]
+	    if ((flag == ALWAYSYES) || (flag == ALWAYSNO))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    if (flag == NO) {
+	        if (y1 < fluxratio) {
+		    Memi[f+i] = YES
+		    call gseti (gp, G_PMLTYPE, 0)
+		    call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+		    call gseti (gp, G_PMLTYPE, 1)
+		    call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+		}
+	    } else {
+	        if (y1 >= fluxratio) {
+		    Memi[f+i] = NO
+		    call gseti (gp, G_PMLTYPE, 0)
+		    call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+		    call gseti (gp, G_PMLTYPE, 1)
+		    call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+		}
+	    }
+	}
+end
+
+
+# CR_PLOT -- Make log plot
+
+procedure cr_plot (cr, im, fluxratio)
+
+pointer	cr			# Cosmic ray list
+pointer	im			# Image pointer
+real	fluxratio		# Flux ratio threshold
+
+int	fd, open(), errcode()
+pointer	sp, fname, gp, gt, gopen(), gt_init()
+errchk	gopen
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	# Open the plotfile.
+	call clgstr ("plotfile", Memc[fname], SZ_FNAME)
+	iferr (fd = open (Memc[fname], APPEND, BINARY_FILE)) {
+	    if (errcode() != SYS_FNOFNAME)
+		call erract (EA_WARN)
+	    return
+	}
+
+	# Set up the graphics.
+	gp = gopen ("stdplot", NEW_FILE, fd)
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_sets (gt, GTXTRAN, "log")
+	call gt_setr (gt, GTXMIN, 10.)
+	call gt_setr (gt, GTYMIN, 0.)
+	call gt_sets (gt, GTTITLE, "Parameters of cosmic rays candidates")
+	call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+	call gt_sets (gt, GTXLABEL, "Flux")
+	call gt_sets (gt, GTYLABEL, "Flux Ratio")
+
+	call cr_graph (gp, gt, cr, fluxratio, 'r')
+
+	call gt_free (gt)
+	call gclose (gp)
+	call close (fd)
+	call sfree (sp)
+end
+
+
+# CR_SHOW -- Select data to show.
+# This returns pointers to the data.  Note the pointers are salloc from
+# the last smark which is done by the calling program.
+
+procedure cr_show (show, cr, x, y, w, flag, index, ncr)
+
+int	show		#I Data to show (0=all, 1=train)
+pointer	cr		#I CR data
+pointer	x		#O Fluxes
+pointer	y		#O Ratios
+pointer	w		#O Weights
+pointer	flag		#O Flags
+pointer	index		#O Index into CR data (if not null)
+int	ncr		#O Number of selected data points
+
+int	i
+
+begin
+	switch (show) {
+	case 0:
+	    ncr = CR_NCR(cr)
+	    x = CR_FLUX(cr) - 1
+	    y = CR_RATIO(cr) - 1
+	    w = CR_WT(cr) - 1
+	    flag = CR_FLAG(cr) - 1
+	    index = NULL
+	case 1:
+	    ncr = CR_NCR(cr)
+	    call salloc (x, ncr, TY_REAL)
+	    call salloc (y, ncr, TY_REAL)
+	    call salloc (w, ncr, TY_REAL)
+	    call salloc (flag, ncr, TY_INT)
+	    call salloc (index, ncr, TY_INT)
+
+	    ncr = 0
+	    x = x - 1
+	    y = y - 1
+	    w = w - 1
+	    flag = flag - 1
+	    index = index - 1
+
+	    do i = 1, CR_NCR(cr) {
+		if (Memr[CR_WT(cr)+i-1] == 0.)
+		    next
+		ncr = ncr + 1
+		Memr[x+ncr] = Memr[CR_FLUX(cr)+i-1]
+		Memr[y+ncr] = Memr[CR_RATIO(cr)+i-1]
+		Memr[w+ncr] = Memr[CR_WT(cr)+i-1]
+		Memi[flag+ncr] = Memi[CR_FLAG(cr)+i-1]
+		Memi[index+ncr] = i
+	    }
+	}
+end
diff --git a/noao/imred/ccdred/src/cosmic/crfind.x b/noao/imred/ccdred/src/cosmic/crfind.x
new file mode 100644
index 00000000..58850940
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/crfind.x
@@ -0,0 +1,305 @@
+include	<math/gsurfit.h>
+
+# CR_FIND -- Find cosmic ray candidates.
+# This procedure is an interface to special procedures specific to a given
+# window size.
+
+procedure cr_find (cr, threshold, data, nc, nl, col, line,
+    sf1, sf2, x, y, z, w)
+
+pointer	cr						# Cosmic ray list
+real	threshold					# Detection threshold
+pointer	data[ARB]					# Data lines
+int	nc						# Number of columns
+int	nl						# Number of lines
+int	col						# First column
+int	line						# Center line
+pointer	sf1, sf2					# Surface fitting
+real	x[ARB], y[ARB], z[ARB], w[ARB]			# Surface arrays
+
+pointer	a, b, c, d, e, f, g
+
+begin
+	switch (nl) {
+	case 5:
+	    a = data[1]
+	    b = data[2]
+	    c = data[3]
+	    d = data[4]
+	    e = data[5]
+	    call cr_find5 (cr, threshold,  col, line, Memr[a], Memr[b],
+		Memr[c], Memr[d], Memr[e], nc, sf1, sf2, x, y, z, w)
+	case 7:
+	    a = data[1]
+	    b = data[2]
+	    c = data[3]
+	    d = data[4]
+	    e = data[5]
+	    f = data[6]
+	    g = data[7]
+	    call cr_find7 (cr, threshold,  col, line, Memr[a], Memr[b],
+		Memr[c], Memr[d], Memr[e], Memr[f], Memr[g], nc,
+		sf1, sf2, x, y, z, w)
+	}
+end
+
+
+# CR_FIND7 -- Find cosmic rays candidates in 7x7 window.
+# This routine finds cosmic rays candidates with the following algorithm.
+#    1. If the pixel is not a local maximum relative to it's 48 neighbors
+#	go on to the next pixel.
+#    2. Identify the next strongest pixel in the 7x7 region.
+#	This suspect pixel is excluded in the following.
+#    2. Compute the flux of the 7x7 region excluding the cosmic ray
+#	candidate and the suspect pixel.
+#    3. The candidate must exceed the average flux per pixel by a specified
+#	threshold.  If not go on to the next pixel.
+#    4. Fit a plane to the border pixels (excluding the suspect pixel).
+#    5. Subtract the background defined by the plane.
+#    6. Determine a replacement value as the average of the four adjacent
+#	pixels (excluding the suspect pixels).
+#    7. Add the pixel to the cosmic ray candidate list.
+
+procedure cr_find7 (cr, threshold, col, line, a, b, c, d, e, f, g, n,
+	sf1, sf2, x, y, z, w)
+
+pointer	cr						# Cosmic ray list
+real	threshold					# Detection threshold
+int	col						# First column
+int	line						# Line
+real	a[ARB], b[ARB], c[ARB], d[ARB]			# Image lines
+real	e[ARB], f[ARB], g[ARB]				# Image lines
+int	n						# Number of columns
+pointer	sf1, sf2					# Surface fitting
+real	x[49], y[49], z[49], w[49]			# Surface arrays
+
+real	bkgd[49] 
+int	i1, i2, i3, i4, i5, i6, i7, j, j1, j2
+real	p, flux, replace, asumr()
+pointer	sf
+
+begin
+	for (i4=4; i4<=n-3; i4=i4+1) {
+	    # Must be local maxima.
+	    p = d[i4]
+	    if (p<a[i4]||p<b[i4]||p<c[i4]||p<e[i4]||p<f[i4]||p<g[i4])
+		next
+	    i1 = i4 - 3
+	    if (p<a[i1]||p<b[i1]||p<c[i1]||p<d[i1]||p<e[i1]||p<f[i1]||p<g[i1])
+		next
+	    i2 = i4 - 2
+	    if (p<a[i2]||p<b[i2]||p<c[i2]||p<d[i2]||p<e[i2]||p<f[i2]||p<g[i2])
+		next
+	    i3 = i4 - 1
+	    if (p<a[i3]||p<b[i3]||p<c[i3]||p<d[i3]||p<e[i3]||p<f[i3]||p<g[i3])
+		next
+	    i5 = i4 + 1
+	    if (p<a[i5]||p<b[i5]||p<c[i5]||p<d[i5]||p<e[i5]||p<f[i5]||p<g[i5])
+		next
+	    i6 = i4 + 2
+	    if (p<a[i6]||p<b[i6]||p<c[i6]||p<d[i6]||p<e[i6]||p<f[i6]||p<g[i6])
+		next
+	    i7 = i4 + 3
+	    if (p<a[i7]||p<b[i7]||p<c[i7]||p<d[i7]||p<e[i7]||p<f[i7]||p<g[i7])
+		next
+
+	    # Convert to a single array in surface fitting order. 
+	    call amovr (a[i1], z[1], 7)
+	    z[8] = b[i7]; z[9] = c[i7]; z[10] = d[i7]; z[11] = e[i7]
+	    z[12] = f[i7]; z[13] = g[i7]; z[14] = g[i6]; z[15] = g[i5]
+	    z[16] = f[i4]; z[17] = g[i3]; z[18] = g[i2]; z[19] = g[i1]
+	    z[20] = f[i1]; z[21] = e[i1]; z[22] = d[i1]; z[23] = c[i1]
+	    z[24] = b[i1]
+	    call amovr (b[i2], z[25], 5)
+	    call amovr (c[i2], z[30], 5)
+	    call amovr (d[i2], z[35], 5)
+	    call amovr (e[i2], z[40], 5)
+	    call amovr (f[i2], z[45], 5)
+
+	    # Find the highest point excluding the center.
+	    j1 = 37; j2 = 1
+	    do j = 2, 49 {
+		if (j == j1)
+		    next
+		if (z[j] > z[j2])
+		    j2 = j 
+	    }
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 49) - z[j1] - z[j2]) / 47
+
+	    # Pixel must be exceed specified threshold.
+	    if (p < flux + threshold)
+		next
+
+	    # Fit and subtract the background.
+	    if (j2 < 25) {
+	        w[j2] = 0
+		sf = sf2
+	        call gsfit (sf, x, y, z, w, 24, WTS_USER, j)
+		w[j2] = 1
+	    } else {
+		sf = sf1
+	        call gsrefit (sf, x, y, z, w, j)
+	    }
+
+	    call gsvector (sf, x, y, bkgd, 49)
+	    call asubr (z, bkgd, z, 49)
+	    p = z[j1]
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 49) - z[j1] - z[j2]) / 47
+
+	    # Determine replacement value from four nearest neighbors again
+	    # excluding the most deviant pixels.
+	    replace = 0
+	    j = 0
+	    if (j2 != 32) {
+		replace = replace + c[i4]
+		j = j + 1
+	    }
+	    if (j2 != 36) {
+		replace = replace + d[i3]
+		j = j + 1
+	    }
+	    if (j2 != 38) {
+		replace = replace + d[i5]
+		j = j + 1
+	    }
+	    if (j2 != 42) {
+		replace = replace + e[i4]
+		j = j + 1
+	    }
+	    replace = replace / j
+
+	    # Add pixel to cosmic ray list.
+	    flux = 100. * flux
+	    call cr_add (cr, col+i4-1, line, flux, flux/p, 0., replace, 0)
+	    i4 = i7
+	}
+end
+
+
+# CR_FIND5 -- Find cosmic rays candidates in 5x5 window.
+# This routine finds cosmic rays candidates with the following algorithm.
+#    1. If the pixel is not a local maximum relative to it's 24 neighbors
+#	go on to the next pixel.
+#    2. Identify the next strongest pixel in the 5x5 region.
+#	This suspect pixel is excluded in the following.
+#    2. Compute the flux of the 5x5 region excluding the cosmic ray
+#	candidate and the suspect pixel.
+#    3. The candidate must exceed the average flux per pixel by a specified
+#	threshold.  If not go on to the next pixel.
+#    4. Fit a plane to the border pixels (excluding the suspect pixel).
+#    5. Subtract the background defined by the plane.
+#    6. Determine a replacement value as the average of the four adjacent
+#	pixels (excluding the suspect pixels).
+#    7. Add the pixel to the cosmic ray candidate list.
+
+procedure cr_find5 (cr, threshold, col, line, a, b, c, d, e, n,
+	sf1, sf2, x, y, z, w)
+
+pointer	cr						# Cosmic ray list
+real	threshold					# Detection threshold
+int	col						# First column
+int	line						# Line
+real	a[ARB], b[ARB], c[ARB], d[ARB], e[ARB]		# Image lines
+int	n						# Number of columns
+pointer	sf1, sf2					# Surface fitting
+real	x[25], y[25], z[25], w[25]			# Surface arrays
+
+real	bkgd[25] 
+int	i1, i2, i3, i4, i5, j, j1, j2
+real	p, flux, replace, asumr()
+pointer	sf
+
+begin
+	for (i3=3; i3<=n-2; i3=i3+1) {
+	    # Must be local maxima.
+	    p = c[i3]
+	    if (p<a[i3]||p<b[i3]||p<d[i3]||p<e[i3])
+		next
+	    i1 = i3 - 2
+	    if (p<a[i1]||p<b[i1]||p<c[i1]||p<d[i1]||p<e[i1])
+		next
+	    i2 = i3 - 1
+	    if (p<a[i2]||p<b[i2]||p<c[i2]||p<d[i2]||p<e[i2])
+		next
+	    i4 = i3 + 1
+	    if (p<a[i4]||p<b[i4]||p<c[i4]||p<d[i4]||p<e[i4])
+		next
+	    i5 = i3 + 2
+	    if (p<a[i5]||p<b[i5]||p<c[i5]||p<d[i5]||p<e[i5])
+		next
+
+	    # Convert to a single array in surface fitting order. 
+	    call amovr (a[i1], z[1], 5)
+	    z[6] = b[i5]; z[7] = c[i5]; z[8] = d[i5]; z[9] = e[i5]
+	    z[10] = e[i4]; z[11] = e[i3]; z[12] = e[i2]; z[13] = e[i1]
+	    z[14] = d[i1]; z[15] = c[i1]; z[16] = b[i1]
+	    call amovr (b[i2], z[17], 3)
+	    call amovr (c[i2], z[20], 3)
+	    call amovr (d[i2], z[23], 3)
+
+	    # Find the highest point excluding the center.
+	    j1 = 21; j2 = 1
+	    do j = 2, 25 {
+		if (j == j1)
+		    next
+		if (z[j] > z[j2])
+		    j2 = j 
+	    }
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 25) - z[j1] - z[j2]) / 23
+
+	    # Pixel must be exceed specified threshold.
+	    if (p < flux + threshold)
+		next
+
+	    # Fit and subtract the background.
+	    if (j2 < 17) {
+	        w[j2] = 0
+		sf = sf2
+	        call gsfit (sf, x, y, z, w, 16, WTS_USER, j)
+		w[j2] = 1
+	    } else {
+		sf = sf1
+	        call gsrefit (sf, x, y, z, w, j)
+	    }
+
+	    call gsvector (sf, x, y, bkgd, 25)
+	    call asubr (z, bkgd, z, 25)
+	    p = z[j1]
+
+	    # Compute the flux excluding the extreme points.
+	    flux = (asumr (z, 25) - z[j1] - z[j2]) / 23
+
+	    # Determine replacement value from four nearest neighbors again
+	    # excluding the most deviant pixels.
+	    replace = 0
+	    j = 0
+	    if (j2 != 18) {
+		replace = replace + b[i3]
+		j = j + 1
+	    }
+	    if (j2 != 20) {
+		replace = replace + c[i2]
+		j = j + 1
+	    }
+	    if (j2 != 22) {
+		replace = replace + c[i4]
+		j = j + 1
+	    }
+	    if (j2 != 24) {
+		replace = replace + d[i3]
+		j = j + 1
+	    }
+	    replace = replace / j
+
+	    # Add pixel to cosmic ray list.
+	    flux = 100. * flux
+	    call cr_add (cr, col+i3-1, line, flux, flux/p, 0., replace, 0)
+	    i3 = i5
+	}
+end
diff --git a/noao/imred/ccdred/src/cosmic/crlist.h b/noao/imred/ccdred/src/cosmic/crlist.h
new file mode 100644
index 00000000..1ed498a7
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/crlist.h
@@ -0,0 +1,17 @@
+define	CR_ALLOC	100		# Allocation block size
+define	CR_LENSTRUCT	9		# Length of structure
+
+define	CR_NCR		Memi[$1]	# Number of cosmic rays
+define	CR_NALLOC	Memi[$1+1]	# Length of cosmic ray list
+define	CR_COL		Memi[$1+2]	# Pointer to columns
+define	CR_LINE		Memi[$1+3]	# Pointer to lines
+define	CR_FLUX		Memi[$1+4]	# Pointer to fluxes
+define	CR_RATIO	Memi[$1+5]	# Pointer to flux ratios
+define	CR_WT		Memi[$1+6]	# Pointer to training weights
+define	CR_REPLACE	Memi[$1+7]	# Pointer to replacement values
+define	CR_FLAG		Memi[$1+8]	# Pointer to rejection flag
+
+define	ALWAYSNO	3
+define	ALWAYSYES	4
+
+define	CR_RMAX		3.		# Maximum radius for matching
diff --git a/noao/imred/ccdred/src/cosmic/crlist.x b/noao/imred/ccdred/src/cosmic/crlist.x
new file mode 100644
index 00000000..e0a8fd5c
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/crlist.x
@@ -0,0 +1,366 @@
+include	<error.h>
+include	<syserr.h>
+include	<gset.h>
+include	"crlist.h"
+
+define	HELP	"noao$lib/scr/cosmicrays.key"
+define	PROMPT	"cosmic ray options"
+
+# CR_OPEN    -- Open cosmic ray list
+# CR_CLOSE   -- Close cosmic ray list
+# CR_ADD     -- Add a cosmic ray candidate to cosmic ray list.
+# CR_TRAIN   -- Set flux ratio threshold from a training set.
+# CR_FINDTHRESH -- Find flux ratio.
+# CR_WEIGHT  -- Compute the training weight at a particular flux ratio.
+# CR_FLAGS   -- Set cosmic ray reject flags.
+# CR_BADPIX  -- Store cosmic rays in bad pixel list.
+# CR_REPLACE -- Replace cosmic rays in image with replacement values.
+
+# CR_OPEN -- Open cosmic ray list
+
+procedure cr_open (cr)
+
+pointer	cr		# Cosmic ray list pointer
+errchk	malloc
+
+begin
+	call malloc (cr, CR_LENSTRUCT, TY_STRUCT)
+	call malloc (CR_COL(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_LINE(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_FLUX(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_RATIO(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_WT(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_REPLACE(cr), CR_ALLOC, TY_REAL)
+	call malloc (CR_FLAG(cr), CR_ALLOC, TY_INT)
+	CR_NCR(cr) = 0
+	CR_NALLOC(cr) = CR_ALLOC
+end
+
+
+# CR_CLOSE -- Close cosmic ray list
+
+procedure cr_close (cr)
+
+pointer	cr		# Cosmic ray list pointer
+
+begin
+	call mfree (CR_COL(cr), TY_REAL)
+	call mfree (CR_LINE(cr), TY_REAL)
+	call mfree (CR_FLUX(cr), TY_REAL)
+	call mfree (CR_RATIO(cr), TY_REAL)
+	call mfree (CR_WT(cr), TY_REAL)
+	call mfree (CR_REPLACE(cr), TY_REAL)
+	call mfree (CR_FLAG(cr), TY_INT)
+	call mfree (cr, TY_STRUCT)
+end
+
+# CR_ADD -- Add a cosmic ray candidate to cosmic ray list.
+
+procedure cr_add (cr, col, line, flux, ratio, wt, replace, flag)
+
+pointer	cr		# Cosmic ray list pointer
+int	col		# Cofluxn
+int	line		# Line
+real	flux		# Luminosity
+real	ratio		# Ratio
+real	wt		# Weight
+real	replace		# Sky value
+int	flag		# Flag value
+
+int	ncr
+errchk	realloc
+
+begin
+	if (CR_NCR(cr) == CR_NALLOC(cr)) {
+	    CR_NALLOC(cr) = CR_NALLOC(cr) + CR_ALLOC
+	    call realloc (CR_COL(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_LINE(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_FLUX(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_RATIO(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_WT(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_REPLACE(cr), CR_NALLOC(cr), TY_REAL)
+	    call realloc (CR_FLAG(cr), CR_NALLOC(cr), TY_INT)
+	}
+
+	ncr = CR_NCR(cr)
+	CR_NCR(cr) = ncr + 1
+	Memr[CR_COL(cr)+ncr] = col
+	Memr[CR_LINE(cr)+ncr] = line
+	Memr[CR_FLUX(cr)+ncr] = flux
+	Memr[CR_RATIO(cr)+ncr] = ratio
+	Memr[CR_WT(cr)+ncr] = wt
+	Memr[CR_REPLACE(cr)+ncr] = replace
+	Memi[CR_FLAG(cr)+ncr] = flag
+end
+
+
+# CR_TRAIN -- Set flux ratio threshold from a training set.
+
+procedure cr_train (cr, gp, gt, im, fluxratio, fname)
+
+pointer	cr		#I Cosmic ray list
+pointer	gp		#I GIO pointer
+pointer	gt		#I GTOOLS pointer
+pointer	im		#I IMIO pointer
+real	fluxratio	#O Flux ratio threshold
+char	fname[ARB]	#I Save file name
+
+char	cmd[10]
+bool	gflag
+real	x, y, y1, y2, w, r, rmin
+int	i, j, n, f, ncr, wcs, key, fd, clgcur(), open(), errcode()
+pointer	col, line, ratio, flux, wt, flag
+
+begin
+	# Open save file
+	iferr (fd = open (fname, APPEND, TEXT_FILE)) {
+	    if (errcode() != SYS_FNOFNAME)
+		call erract (EA_WARN)
+	    fd = 0
+	}
+
+	ncr = CR_NCR(cr)
+	col = CR_COL(cr) - 1
+	line = CR_LINE(cr) - 1
+	flux = CR_FLUX(cr) - 1
+	ratio = CR_RATIO(cr) - 1
+	wt = CR_WT(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	gflag = false
+	n = 0
+	while (clgcur ("objects", x, y, wcs, key, cmd, 10) != EOF) {
+	    switch (key) {
+	    case '?':
+		call gpagefile (gp, HELP, PROMPT)
+		next
+	    case 'q':
+		break
+	    case 's':
+		w = 1
+		f = ALWAYSNO
+	    case 'c':
+		w = -1
+		f = ALWAYSYES
+	    case 'g':
+		if (gflag)
+		    call cr_examine (cr, gp, gt, im, fluxratio, 'z')
+		else {
+		    if (n > 1)
+			call cr_findthresh (cr, fluxratio)
+		    call cr_flags (cr, fluxratio)
+		    call cr_examine (cr, gp, gt, im, fluxratio, 'r')
+		    gflag = true
+		}
+		next
+	    default:
+		next
+	    }
+
+	    y1 = y - CR_RMAX
+	    y2 = y + CR_RMAX
+	    for (i=10; i<ncr && y1>Memr[line+i]; i=i+10)
+		;
+	    j = i - 9
+	    rmin = (Memr[col+j] - x) ** 2 + (Memr[line+j] - y) ** 2
+	    for (i=j+1; i<ncr && y2>Memr[line+i]; i=i+1) {
+		r = (Memr[col+i] - x) ** 2 + (Memr[line+i] - y) ** 2
+		if (r < rmin) {
+		    rmin = r
+		    j = i
+		}
+	    }
+	    if (sqrt (rmin) > CR_RMAX)
+		next
+
+	    Memr[wt+j] = w
+	    Memi[flag+j] = f
+	    n = n + 1
+
+	    if (gflag) {
+		if (n > 1) {
+		    call cr_findthresh (cr, r)
+		    call cr_update (gp, r, cr, fluxratio, 0)
+		}
+		call gmark (gp, Memr[flux+j], Memr[ratio+j], GM_BOX, 2., 2.)
+	    }
+	    if (fd > 0) {
+		call fprintf (fd, "%g %g %d %c\n")
+		    call pargr (x)
+		    call pargr (y)
+		    call pargi (wcs)
+		    call pargi (key)
+	    }
+	}
+
+	if (fd > 0)
+	    call close (fd)
+end
+
+
+# CR_FINDTHRESH -- Find flux ratio.
+
+procedure cr_findthresh (cr, fluxratio)
+
+pointer	cr		#I Cosmic ray list
+real	fluxratio	#O Flux ratio threshold
+
+real	w, r, rmin, cr_weight()
+int	i, ncr
+pointer	ratio, wt
+
+begin
+	ncr = CR_NCR(cr)
+	ratio = CR_RATIO(cr) - 1
+	wt = CR_WT(cr) - 1
+
+	fluxratio = Memr[ratio+1]
+	rmin = cr_weight (fluxratio, Memr[ratio+1], Memr[wt+1], ncr)
+	do i = 2, ncr {
+	    if (Memr[wt+i] == 0.)
+		next
+	    r = Memr[ratio+i]
+	    w = cr_weight (r, Memr[ratio+1], Memr[wt+1], ncr)
+	    if (w <= rmin) {
+		if (w == rmin)
+		    fluxratio = min (fluxratio, r)
+		else {
+		    rmin = w
+		    fluxratio = r
+		}
+	    }
+	}
+end
+
+
+# CR_WEIGHT -- Compute the training weight at a particular flux ratio.
+
+real procedure cr_weight (fluxratio, ratio, wts, ncr)
+
+real	fluxratio		#I Flux ratio
+real	ratio[ARB]		#I Ratio Values
+real	wts[ARB]		#I Weights
+int	ncr			#I Number of ratio values
+real	wt			#O Sum of weights
+
+int	i
+
+begin
+	wt = 0.
+	do i = 1, ncr {
+	    if (ratio[i] > fluxratio) {
+		if (wts[i] < 0.)
+		    wt = wt - wts[i]
+	    } else {
+		if (wts[i] > 0.)
+		    wt = wt + wts[i]
+	    }
+	}
+	return (wt)
+end
+
+
+# CR_FLAGS -- Set cosmic ray reject flags.
+
+procedure cr_flags (cr, fluxratio)
+
+pointer	cr			# Cosmic ray candidate list
+real	fluxratio		# Rejection limits
+
+int	i, ncr
+pointer	ratio, flag
+
+begin
+	ncr = CR_NCR(cr)
+	ratio = CR_RATIO(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == ALWAYSYES) || (Memi[flag+i] == ALWAYSNO))
+		next
+	    if (Memr[ratio+i] > fluxratio)
+		Memi[flag+i] = NO
+	    else
+		Memi[flag+i] = YES
+	}
+end
+
+
+# CR_BADPIX -- Store cosmic rays in bad pixel list.
+# This is currently a temporary measure until a real bad pixel list is
+# implemented.
+
+procedure cr_badpix (cr, fname)
+
+pointer	cr		# Cosmic ray list
+char	fname[ARB]	# Bad pixel file name
+
+int	i, ncr, c, l, f, fd, open(), errcode()
+pointer	col, line, ratio, flux, flag
+errchk	open
+
+begin
+	# Open bad pixel file
+	iferr (fd = open (fname, APPEND, TEXT_FILE)) {
+	    if (errcode() != SYS_FNOFNAME)
+		call erract (EA_WARN)
+	    return
+	}
+
+	ncr = CR_NCR(cr)
+	col = CR_COL(cr) - 1
+	line = CR_LINE(cr) - 1
+	flux = CR_FLUX(cr) - 1
+	ratio = CR_RATIO(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	do i = 1, ncr {
+	    f = Memi[flag+i]
+	    if ((f == NO) || (f == ALWAYSNO))
+		next
+
+	    c = Memr[col+i]
+	    l = Memr[line+i]
+	    call fprintf (fd, "%d %d\n")
+	        call pargi (c)
+	        call pargi (l)
+	}
+	call close (fd)
+end
+
+
+# CR_REPLACE -- Replace cosmic rays in image with replacement values.
+
+procedure cr_replace (cr, offset, im, nreplaced)
+
+pointer	cr		# Cosmic ray list
+int	offset		# Offset in list
+pointer	im		# IMIO pointer of output image
+int	nreplaced	# Number replaced (for log)
+
+int	i, ncr, c, l, f
+real	r
+pointer	col, line, replace, flag, imps2r()
+
+begin
+	ncr = CR_NCR(cr)
+	if (ncr <= offset)
+	    return
+
+	col = CR_COL(cr) - 1
+	line = CR_LINE(cr) - 1
+	replace = CR_REPLACE(cr) - 1
+	flag = CR_FLAG(cr) - 1
+
+	do i = offset+1, ncr {
+	    f = Memi[flag+i]
+	    if ((f == NO) || (f == ALWAYSNO))
+		next
+
+	    c = Memr[col+i]
+	    l = Memr[line+i]
+	    r = Memr[replace+i]
+	    Memr[imps2r (im, c, c, l, l)] = r
+	    nreplaced = nreplaced + 1
+	}
+end
diff --git a/noao/imred/ccdred/src/cosmic/crsurface.x b/noao/imred/ccdred/src/cosmic/crsurface.x
new file mode 100644
index 00000000..32645ff4
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/crsurface.x
@@ -0,0 +1,46 @@
+define	DUMMY	6
+
+# CR_SURFACE -- Draw a perspective view of a surface.  The altitude
+# and azimuth of the viewing angle are variable.
+
+procedure cr_surface(gp, data, ncols, nlines, angh, angv)
+
+pointer	gp			# GIO pointer
+real	data[ncols,nlines]	# Surface data to be plotted
+int	ncols, nlines		# Dimensions of surface
+real	angh, angv		# Orientation of surface (degrees)
+
+int	wkid
+pointer	sp, work
+
+int	first
+real	vpx1, vpx2, vpy1, vpy2
+common  /frstfg/ first
+common  /noaovp/ vpx1, vpx2, vpy1, vpy2
+
+begin
+	call smark (sp)
+	call salloc (work, 2 * (2 * ncols * nlines + ncols + nlines), TY_REAL)
+
+	# Initialize surface common blocks
+	first = 1
+	call srfabd()
+
+	# Define viewport.
+	call ggview (gp, vpx1, vpx2, vpy1, vpy2) 
+
+	# Link GKS to GIO
+	wkid = 1
+	call gopks (STDERR)
+	call gopwk (wkid, DUMMY, gp)
+	call gacwk (wkid)
+
+	call ezsrfc (data, ncols, nlines, angh, angv, Memr[work])
+
+	call gdawk (wkid)
+	# We don't want to close the GIO pointer.
+	#call gclwk (wkid)
+	call gclks ()
+
+	call sfree (sp)
+end
diff --git a/noao/imred/ccdred/src/cosmic/mkpkg b/noao/imred/ccdred/src/cosmic/mkpkg
new file mode 100644
index 00000000..d63d9c2c
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/mkpkg
@@ -0,0 +1,16 @@
+# COSMIC RAY CLEANING
+
+$checkout libpkg.a ../..
+$update   libpkg.a
+$checkin  libpkg.a ../..
+$exit
+
+libpkg.a:
+	crexamine.x	crlist.h <error.h> <gset.h> <mach.h> <pkg/gtools.h>\
+			<imhdr.h> <syserr.h>
+	crfind.x	<math/gsurfit.h>
+	crlist.x	crlist.h <error.h> <gset.h> <syserr.h>
+	crsurface.x	
+	t_cosmicrays.x	crlist.h <error.h> <gset.h> <math/gsurfit.h>\
+			<pkg/gtools.h> <imhdr.h> <imset.h>
+	;
diff --git a/noao/imred/ccdred/src/cosmic/t_cosmicrays.x b/noao/imred/ccdred/src/cosmic/t_cosmicrays.x
new file mode 100644
index 00000000..8640b639
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/t_cosmicrays.x
@@ -0,0 +1,348 @@
+include	<error.h>
+include <imhdr.h>
+include <imset.h>
+include	<math/gsurfit.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	"crlist.h"
+
+# T_COSMICRAYS -- Detect and remove cosmic rays in images.
+# A list of images is examined for cosmic rays which are then replaced
+# by values from neighboring pixels.  The output image may be the same
+# as the input image.  This is the top level procedure which manages
+# the input and output image data.  The actual algorithm for detecting
+# cosmic rays is in CR_FIND.
+
+procedure t_cosmicrays ()
+
+int	list1			# List of input images to be cleaned
+int	list2			# List of output images
+int	list3			# List of output bad pixel files
+real	threshold		# Detection threshold
+real	fluxratio		# Luminosity boundary for stars
+int	npasses			# Number of cleaning passes
+int	szwin			# Size of detection window
+bool	train			# Use training objects?
+pointer	savefile		# Save file for training objects
+bool	interactive		# Examine cosmic ray parameters?
+char	ans			# Answer to interactive query
+
+int	nc, nl, c, c1, c2, l, l1, l2, szhwin, szwin2
+int	i, j, k, m, ncr, ncrlast, nreplaced, flag
+pointer	sp, input, output, badpix, str, gp, gt, im, in, out
+pointer	x, y, z, w, sf1, sf2, cr, data, ptr
+
+bool	clgetb(), ccdflag(), streq(), strne()
+char	clgetc()
+int	imtopenp(), imtlen(), imtgetim(), clpopnu(), clgfil(), clgeti()
+real	clgetr()
+pointer	immap(), impl2r(), imgs2r(), gopen(), gt_init()
+errchk	immap, impl2r, imgs2r
+errchk	cr_find, cr_examine, cr_replace, cr_plot, cr_badpix
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (badpix, SZ_FNAME, TY_CHAR)
+	call salloc (savefile, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the task parameters.  Check that the number of output images
+	# is either zero, in which case the cosmic rays will be removed
+	# in place, or equal to the number of input images.
+
+	list1 = imtopenp ("input")
+	list2 = imtopenp ("output")
+	i = imtlen (list1)
+	j = imtlen (list2)
+	if (j > 0 && j != i)
+	    call error (0, "Input and output image lists do not match")
+
+	list3 = clpopnu ("badpix")
+	threshold = clgetr ("threshold")
+	fluxratio = clgetr ("fluxratio")
+	npasses = clgeti ("npasses")
+	szwin = clgeti ("window")
+	train = clgetb ("train")
+	call clgstr ("savefile", Memc[savefile], SZ_FNAME)
+	interactive = clgetb ("interactive")
+	call clpstr ("answer", "yes")
+	ans = 'y'
+
+	# Set up the graphics.
+	call clgstr ("graphics", Memc[str], SZ_LINE)
+	if (interactive) {
+	    gp = gopen (Memc[str], NEW_FILE+AW_DEFER, STDGRAPH)
+	    gt = gt_init()
+	    call gt_sets (gt, GTTYPE, "mark")
+	    call gt_sets (gt, GTXTRAN, "log")
+	    call gt_setr (gt, GTXMIN, 10.)
+	    call gt_setr (gt, GTYMIN, 0.)
+	    call gt_sets (gt, GTTITLE, "Parameters of cosmic rays candidates")
+	    call gt_sets (gt, GTXLABEL, "Flux")
+	    call gt_sets (gt, GTYLABEL, "Flux Ratio")
+	}
+
+	# Use image header translation file.
+	call clgstr ("instrument", Memc[input], SZ_FNAME)
+	call hdmopen (Memc[input])
+
+	# Set up surface fitting.  The background points are placed together
+	# at the beginning of the arrays.  There are two surface pointers,
+	# one for using the fast refit if there are no points excluded and
+	# one for doing a full fit with points excluded.
+
+	szhwin = szwin / 2
+	szwin2 = szwin * szwin
+	call salloc (data, szwin, TY_INT)
+	call salloc (x, szwin2, TY_REAL)
+	call salloc (y, szwin2, TY_REAL)
+	call salloc (z, szwin2, TY_REAL)
+	call salloc (w, szwin2, TY_REAL)
+
+	k = 0
+	do i = 1, szwin {
+	    Memr[x+k] = i
+	    Memr[y+k] = 1
+	    k = k + 1
+	}
+	do i = 2, szwin {
+	    Memr[x+k] = szwin
+	    Memr[y+k] = i
+	    k = k + 1
+	}
+	do i = szwin-1, 1, -1 {
+	    Memr[x+k] = i
+	    Memr[y+k] = szwin
+	    k = k + 1
+	}
+	do i = szwin-1, 2, -1 {
+	    Memr[x+k] = 1
+	    Memr[y+k] = i
+	    k = k + 1
+	}
+	do i = 2, szwin-1 {
+	    do j = 2, szwin-1 {
+	        Memr[x+k] = j
+	        Memr[y+k] = i
+	        k = k + 1
+	    }
+	}
+	call aclrr (Memr[z], szwin2)
+	call amovkr (1., Memr[w], 4*szwin-4)
+	call gsinit (sf1, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin),
+	    1., real(szwin))
+	call gsinit (sf2, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin),
+	    1., real(szwin))
+	call gsfit (sf1, Memr[x], Memr[y], Memr[z], Memr[w], 4*szwin-4,
+	    WTS_USER, j)
+
+	# Process each input image.  Either work in place or create a
+	# new output image.  If an error mapping the images occurs
+	# issue a warning and go on to the next input image.
+
+	while (imtgetim (list1, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (list2, Memc[output], SZ_FNAME) == EOF)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+	    if (clgfil (list3, Memc[badpix], SZ_FNAME) == EOF)
+		Memc[badpix] = EOS
+
+	    iferr {
+		in = NULL
+		out = NULL
+		cr = NULL
+
+		# Map the input image and check for image type and
+		# previous correction flag.  If the output image is
+		# the same as the input image work in place.
+		# Initialize IMIO to use a scrolling buffer of lines.
+
+		call set_input (Memc[input], im, i)
+		if (im == NULL)
+		    call error (1, "Skipping input image")
+
+		if (ccdflag (im, "crcor")) {
+		    call eprintf ("WARNING: %s previously corrected\n")
+			call pargstr (Memc[input])
+		    #call imunmap (im)
+		    #next
+		}
+
+	        if (streq (Memc[input], Memc[output])) {
+		    call imunmap (im)
+	            im = immap (Memc[input], READ_WRITE, 0)
+		}
+		in = im
+
+	        nc = IM_LEN(in,1)
+	        nl = IM_LEN(in,2)
+	        if ((nl < szwin) || (nc < szwin))
+	            call error (0, "Image size is too small")
+	        call imseti (in, IM_NBUFS, szwin)
+	        call imseti (in, IM_TYBNDRY, BT_NEAREST)
+	        call imseti (in, IM_NBNDRYPIX, szhwin)
+
+		# Open the output image if needed.
+	        if (strne (Memc[input], Memc[output]))
+		    im = immap (Memc[output], NEW_COPY, in)
+		out = im
+
+		# Open a cosmic ray list structure.
+	        call cr_open (cr)
+		ncrlast = 0
+		nreplaced = 0
+
+		# Now proceed through the image line by line, scrolling
+		# the line buffers at each step.  If creating a new image
+		# also write out each line as it is read.  A procedure is
+		# called to find the cosmic ray candidates in the line
+		# and add them to the list maintained by CRLIST.
+		# Note that cosmic rays are not replaced at this point
+		# in order to allow the user to modify the criteria for
+		# a cosmic ray and review the results.
+
+	        c1 = 1-szhwin
+	        c2 = nc+szhwin
+		do i = 1, szwin-1 
+		    Memi[data+i] =
+			imgs2r (in, c1, c2, i-szhwin, i-szhwin)
+
+	        do l = 1, nl {
+		    do i = 1, szwin-1
+			Memi[data+i-1] = Memi[data+i]
+		    Memi[data+szwin-1] =
+			imgs2r (in, c1, c2, l+szhwin, l+szhwin)
+		    if (out != in)
+		        call amovr (Memr[Memi[data+szhwin]+szhwin],
+			    Memr[impl2r(out,l)], nc)
+
+	            call cr_find (cr, threshold, Memi[data],
+		        c2-c1+1, szwin, c1, l,
+		        sf1, sf2, Memr[x], Memr[y], Memr[z], Memr[w])
+	        }
+		if (interactive && train) {
+		    call cr_train (cr, gp, gt, in, fluxratio, Memc[savefile])
+		    train = false
+		}
+	        call cr_flags (cr, fluxratio)
+
+		# If desired examine the cosmic ray list interactively.
+	        if (interactive && ans != 'N') {
+		    if (ans != 'Y') {
+		        call eprintf ("%s - ")
+		            call pargstr (Memc[input])
+		        call flush (STDERR)
+		        ans = clgetc ("answer")
+		    }
+		    if ((ans == 'Y') || (ans == 'y'))
+	                call cr_examine (cr, gp, gt, in, fluxratio, 'r')
+	        }
+
+		# Now replace the selected cosmic rays in the output image.
+
+		call imflush (out)
+		call imseti (out, IM_ADVICE, RANDOM)
+	        call cr_replace (cr, ncrlast, out, nreplaced)
+
+		# Do additional passes through the data.  We work in place
+		# in the output image.  Note that we only have to look in
+		# the vicinity of replaced cosmic rays for secondary
+		# events since we've already looked at every pixel once.
+		# Instead of scrolling through the image we will extract
+		# subrasters around each replaced cosmic ray.  However,
+		# we use pointers into the subraster to maintain the same
+		# format expected by CR_FIND.
+
+	        if (npasses > 1) {
+	            if (out != in)
+	                call imunmap (out)
+	            call imunmap (in)
+		    im = immap (Memc[output], READ_WRITE, 0)
+		    in = im
+		    out = im
+	            call imseti (in, IM_TYBNDRY, BT_NEAREST)
+	            call imseti (in, IM_NBNDRYPIX, szhwin)
+
+	            for (i=2; i<=npasses; i=i+1) {
+			# Loop through each cosmic ray in the previous pass.
+		        ncr = CR_NCR(cr)
+		        do j = ncrlast+1, ncr {
+			    flag = Memi[CR_FLAG(cr)+j-1]
+			    if (flag==NO || flag==ALWAYSNO)
+			        next
+			    c = Memr[CR_COL(cr)+j-1]
+			    l = Memr[CR_LINE(cr)+j-1]
+			    c1 = max (1-szhwin, c - (szwin-1))
+			    c2 = min (nc+szhwin, c + (szwin-1))
+			    k = c2 - c1 + 1
+			    l1 = max (1-szhwin, l - (szwin-1))
+			    l2 = min (nl+szhwin, l + (szwin-1))
+
+			    # Set the line pointers off an image section
+			    # centered on a previously replaced cosmic ray.
+
+			    ptr = imgs2r (in, c1, c2, l1, l2) - k
+
+			    l1 = max (1, l - szhwin)
+			    l2 = min (nl, l + szhwin)
+			    do l = l1, l2 {
+				do m = 1, szwin
+				    Memi[data+m-1] = ptr + m * k
+				ptr = ptr + k
+
+	                        call cr_find ( cr, threshold, Memi[data],
+				    k, szwin, c1, l, sf1, sf2,
+				    Memr[x], Memr[y], Memr[z], Memr[w])
+			    }
+	                }
+		        call cr_flags (cr, fluxratio)
+
+			# Replace any new cosmic rays found.
+	                call cr_replace (cr, ncr, in, nreplaced)
+			ncrlast = ncr
+		    }
+	        }
+
+		# Output header log, log, plot, and bad pixels.
+		call sprintf (Memc[str], SZ_LINE,
+		    "Threshold=%5.1f, fluxratio=%6.2f, removed=%d")
+		    call pargr (threshold)
+		    call pargr (fluxratio)
+		    call pargi (nreplaced)
+		call timelog (Memc[str], SZ_LINE)
+		call ccdlog (out, Memc[str])
+		call hdmpstr (out, "crcor", Memc[str])
+
+		call cr_plot (cr, in, fluxratio)
+		call cr_badpix (cr, Memc[badpix])
+
+	        call cr_close (cr)
+	        if (out != in)
+	            call imunmap (out)
+	        call imunmap (in)
+	    } then {
+		# In case of error clean up and go on to the next image.
+		if (in != NULL) {
+		    if (out != NULL && out != in)
+			call imunmap (out)
+		    call imunmap (in)
+		}
+		if (cr != NULL)
+		    call cr_close (cr)
+		call erract (EA_WARN)
+	    }
+	}
+
+	if (interactive) {
+	    call gt_free (gt)
+	    call gclose (gp)
+	}
+	call imtclose (list1)
+	call imtclose (list2)
+	call clpcls (list3)
+	call hdmclose ()
+	call gsfree (sf1)
+	call gsfree (sf2)
+	call sfree (sp)
+end