Repatch (from linux) of OSX IRAF

5 files changed, 1976 insertions, 0 deletions

diff --git a/noao/imred/quadred/src/ccdproc/generic/ccdred.h b/noao/imred/quadred/src/ccdproc/generic/ccdred.h
new file mode 100644
index 00000000..ef41f592
--- /dev/null
+++ b/noao/imred/quadred/src/ccdproc/generic/ccdred.h
@@ -0,0 +1,155 @@
+# CCDRED Data Structures and Definitions
+
+# The CCD structure:  This structure is used to communicate processing
+# parameters between the package procedures.  It contains pointers to
+# data, calibration image IMIO pointers, scaling parameters, and the
+# correction flags.  The corrections flags indicate which processing
+# operations are to be performed.  The subsection parameters do not
+# include a step size.  A step size is assumed.  If arbitrary subsampling
+# is desired this would be the next generalization.
+
+define	LEN_CCD		75		# Length of CCD structure
+
+# CCD data coordinates
+define	CCD_C1		Memi[$1]	# CCD starting column
+define	CCD_C2		Memi[$1+1]	# CCD ending column
+define	CCD_L1		Memi[$1+2]	# CCD starting line
+define	CCD_L2		Memi[$1+3]	# CCD ending line
+
+# Input data
+define	IN_IM		Memi[$1+4]	# Input image pointer
+define	IN_C1		Memi[$1+5]	# Input data starting column
+define	IN_C2		Memi[$1+6]	# Input data ending column
+define	IN_L1		Memi[$1+7]	# Input data starting line
+define	IN_L2		Memi[$1+8]	# Input data ending line
+define	IN_NSEC		Memi[$1+71]	# Number of input pieces
+define	IN_SEC		Memi[$1+72]	# Pointer to sections (c1,c2,l1,l2)xn
+
+# Output data
+define	OUT_IM		Memi[$1+9]	# Output image pointer
+define	OUT_C1		Memi[$1+10]	# Output data starting column
+define	OUT_C2		Memi[$1+11]	# Output data ending column
+define	OUT_L1		Memi[$1+12]	# Output data starting line
+define	OUT_L2		Memi[$1+13]	# Output data ending line
+define	OUT_SEC		Memi[$1+73]	# Pointer to sections (c1,c2,l1,l2)xn
+
+# Zero level data
+define	ZERO_IM		Memi[$1+14]	# Zero level image pointer
+define	ZERO_C1		Memi[$1+15]	# Zero level data starting column
+define	ZERO_C2		Memi[$1+16]	# Zero level data ending column
+define	ZERO_L1		Memi[$1+17]	# Zero level data starting line
+define	ZERO_L2		Memi[$1+18]	# Zero level data ending line
+
+# Dark count data
+define	DARK_IM		Memi[$1+19]	# Dark count image pointer
+define	DARK_C1		Memi[$1+20]	# Dark count data starting column
+define	DARK_C2		Memi[$1+21]	# Dark count data ending column
+define	DARK_L1		Memi[$1+22]	# Dark count data starting line
+define	DARK_L2		Memi[$1+23]	# Dark count data ending line
+
+# Flat field data
+define	FLAT_IM		Memi[$1+24]	# Flat field image pointer
+define	FLAT_C1		Memi[$1+25]	# Flat field data starting column
+define	FLAT_C2		Memi[$1+26]	# Flat field data ending column
+define	FLAT_L1		Memi[$1+27]	# Flat field data starting line
+define	FLAT_L2		Memi[$1+28]	# Flat field data ending line
+
+# Illumination data
+define	ILLUM_IM	Memi[$1+29]	# Illumination image pointer
+define	ILLUM_C1	Memi[$1+30]	# Illumination data starting column
+define	ILLUM_C2	Memi[$1+31]	# Illumination data ending column
+define	ILLUM_L1	Memi[$1+32]	# Illumination data starting line
+define	ILLUM_L2	Memi[$1+33]	# Illumination data ending line
+
+# Fringe data
+define	FRINGE_IM	Memi[$1+34]	# Fringe image pointer
+define	FRINGE_C1	Memi[$1+35]	# Fringe data starting column
+define	FRINGE_C2	Memi[$1+36]	# Fringe data ending column
+define	FRINGE_L1	Memi[$1+37]	# Fringe data starting line
+define	FRINGE_L2	Memi[$1+38]	# Fringe data ending line
+
+# Trim section
+define	TRIM_C1		Memi[$1+39]	# Trim starting column
+define	TRIM_C2		Memi[$1+40]	# Trim ending column
+define	TRIM_L1		Memi[$1+41]	# Trim starting line
+define	TRIM_L2		Memi[$1+42]	# Trim ending line
+
+# Bias section
+define	BIAS_C1		Memi[$1+43]	# Bias starting column
+define	BIAS_C2		Memi[$1+44]	# Bias ending column
+define	BIAS_L1		Memi[$1+45]	# Bias starting line
+define	BIAS_L2		Memi[$1+46]	# Bias ending line
+define	BIAS_SEC	Memi[$1+74]	# Multiple bias sections
+
+define	READAXIS	Memi[$1+47]	# Read out axis (1=cols, 2=lines)
+define	CALCTYPE	Memi[$1+48]	# Calculation data type
+define	NBADCOLS	Memi[$1+49]	# Number of column interpolation regions
+define	BADCOLS		Memi[$1+50]	# Pointer to col interpolation regions
+define	NBADLINES	Memi[$1+51]	# Number of line interpolation regions
+define	BADLINES	Memi[$1+52]	# Pointer to line interpolation regions
+define	OVERSCAN_VEC	Memi[$1+53]	# Pointer to overscan vector
+define	DARKSCALE	Memr[P2R($1+54)]	# Dark count scale factor
+define	FRINGESCALE	Memr[P2R($1+55)]	# Fringe scale factor
+define	FLATSCALE	Memr[P2R($1+56)]	# Flat field scale factor
+define	ILLUMSCALE	Memr[P2R($1+57)]	# Illumination scale factor
+define	MINREPLACE	Memr[P2R($1+58)]	# Minimum replacement value
+define	MEAN		Memr[P2R($1+59)]	# Mean of output image
+define	COR		Memi[$1+60]	# Overall correction flag
+define	CORS		Memi[$1+61+($2-1)]  # Individual correction flags
+
+# Individual components of input, output, and bias section pieces.
+define	IN_SC1		Memi[IN_SEC($1)+4*$2-4]
+define	IN_SC2		Memi[IN_SEC($1)+4*$2-3]
+define	IN_SL1		Memi[IN_SEC($1)+4*$2-2]
+define	IN_SL2		Memi[IN_SEC($1)+4*$2-1]
+define	OUT_SC1		Memi[OUT_SEC($1)+4*$2-4]
+define	OUT_SC2		Memi[OUT_SEC($1)+4*$2-3]
+define	OUT_SL1		Memi[OUT_SEC($1)+4*$2-2]
+define	OUT_SL2		Memi[OUT_SEC($1)+4*$2-1]
+define	BIAS_SC1	Memi[BIAS_SEC($1)+4*$2-4]
+define	BIAS_SC2	Memi[BIAS_SEC($1)+4*$2-3]
+define	BIAS_SL1	Memi[BIAS_SEC($1)+4*$2-2]
+define	BIAS_SL2	Memi[BIAS_SEC($1)+4*$2-1]
+
+# The correction array contains the following elements with array indices
+# given by the macro definitions.
+
+define	NCORS		10		# Number of corrections
+
+define	FIXPIX		1		# Fix bad pixels
+define	TRIM		2		# Trim image
+define	OVERSCAN	3		# Apply overscan correction
+define	ZEROCOR		4		# Apply zero level correction
+define	DARKCOR		5		# Apply dark count correction
+define	FLATCOR		6		# Apply flat field correction
+define	ILLUMCOR	7		# Apply illumination correction
+define	FRINGECOR	8		# Apply fringe correction
+define	FINDMEAN	9		# Find the mean of the output image
+define	MINREP		10		# Check and replace minimum value
+
+# The following definitions identify the correction values in the correction
+# array.  They are defined in terms of bit fields so that it is possible to
+# add corrections to form unique combination corrections.  Some of
+# these combinations are implemented as compound operations for efficiency.
+
+define	O	001B	# overscan
+define	Z	002B	# zero level
+define	D	004B	# dark count
+define	F	010B	# flat field
+define	I	020B	# Illumination
+define	Q	040B	# Fringe
+
+# The following correction combinations are recognized.
+
+define	ZO	003B	# zero level + overscan
+define	DO	005B	# dark count + overscan
+define	DZ	006B	# dark count + zero level
+define	DZO	007B	# dark count + zero level + overscan
+define	FO	011B	# flat field + overscan
+define	FZ	012B	# flat field + zero level
+define	FZO	013B	# flat field + zero level + overscan
+define	FD	014B	# flat field + dark count
+define	FDO	015B	# flat field + dark count + overscan
+define	FDZ	016B	# flat field + dark count + zero level
+define	FDZO	017B	# flat field + dark count + zero level + overscan
+define	QI	060B	# fringe + illumination
diff --git a/noao/imred/quadred/src/ccdproc/generic/cor.x b/noao/imred/quadred/src/ccdproc/generic/cor.x
new file mode 100644
index 00000000..0dc21310
--- /dev/null
+++ b/noao/imred/quadred/src/ccdproc/generic/cor.x
@@ -0,0 +1,695 @@
+include	"ccdred.h"
+
+
+.help cor Feb87 noao.imred.ccdred
+.nf ----------------------------------------------------------------------------
+cor -- Process CCD image lines
+
+These procedures are the heart of the CCD processing.  They do the desired
+set of processing operations on the image line data as efficiently as
+possible.  They are called by the PROC procedures.  There are four procedures
+one for each readout axis and one for short and real image data.
+Some sets of operations are coded as single compound operations for efficiency.
+To keep the number of combinations managable only the most common
+combinations are coded as compound operations.  The combinations
+consist of any set of line overscan, column overscan, zero level, dark
+count, and flat field and any set of illumination and fringe
+correction.  The corrections are applied in place to the output vector.
+
+The column readout procedure is more complicated in order to handle
+zero level and flat field corrections specified as one dimensional
+readout corrections instead of two dimensional calibration images.
+Column readout format is probably extremely rare and the 1D readout
+corrections are used only for special types of data.
+.ih
+SEE ALSO
+proc, ccdred.h
+.endhelp -----------------------------------------------------------------------
+
+
+# COR1 -- Correct image lines with readout axis 1 (lines).
+
+procedure cor1s (cors, out, overscan, zero, dark, flat, illum,
+	fringe, n, darkscale, flatscale, illumscale, frgscale)
+
+int	cors[ARB]	# Correction flags
+short	out[n]		# Output data
+real	overscan	# Overscan value
+short	zero[n]		# Zero level correction
+short	dark[n]		# Dark count correction
+short	flat[n]		# Flat field correction
+short	illum[n]	# Illumination correction
+short	fringe[n]	# Fringe correction
+int	n		# Number of pixels
+real	darkscale	# Dark count scale factor
+real	flatscale	# Flat field scale factor
+real	illumscale	# Illumination scale factor
+real	frgscale	# Fringe scale factor
+
+int	i, op
+
+begin
+	op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+	switch (op) {
+	case O: # overscan
+	    do i = 1, n
+	        out[i] = out[i] - overscan
+	case Z: # zero level
+	    do i = 1, n
+		out[i] = out[i] - zero[i]
+
+	case ZO: # zero level + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan - zero[i]
+
+	case D: # dark count
+	    do i = 1, n
+		out[i] = out[i] - darkscale * dark[i]
+	case DO: # dark count + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan - darkscale * dark[i]
+	case DZ: # dark count + zero level
+	    do i = 1, n
+		out[i] = out[i] - zero[i] - darkscale * dark[i]
+	case DZO: # dark count + zero level + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan - zero[i] - darkscale * dark[i]
+
+	case F: # flat field
+	    do i = 1, n
+		out[i] = out[i] * flatscale / flat[i]
+	case FO: # flat field + overscan
+	    do i = 1, n
+	        out[i] = (out[i] - overscan) * flatscale / flat[i]
+	case FZ: # flat field + zero level
+	    do i = 1, n
+		out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+	case FZO: # flat field + zero level + overscan
+	    do i = 1, n
+		out[i] = (out[i] - overscan - zero[i]) * flatscale /
+		    flat[i]
+	case FD: # flat field + dark count
+	    do i = 1, n
+		out[i] = (out[i] - darkscale * dark[i]) * flatscale / flat[i]
+	case FDO: # flat field + dark count + overscan
+	    do i = 1, n
+		out[i] = (out[i] - overscan - darkscale * dark[i]) *
+		    flatscale / flat[i]
+	case FDZ: # flat field + dark count + zero level
+	    do i = 1, n
+		out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+		    flatscale / flat[i]
+	case FDZO: # flat field + dark count + zero level + overscan
+	    do i = 1, n
+		out[i] = (out[i] - overscan - zero[i] -
+		    darkscale * dark[i]) * flatscale / flat[i]
+	}
+
+	# Often these operations will not be performed so test for no
+	# correction rather than go through the switch.
+
+	op = cors[ILLUMCOR] + cors[FRINGECOR]
+	if (op != 0) {
+	    switch (op) {
+	    case I: # illumination
+	        do i = 1, n
+		    out[i] = out[i] * illumscale / illum[i]
+	    case Q: # fringe
+	        do i = 1, n
+		    out[i] = out[i] - frgscale * fringe[i]
+	    case QI: # fringe + illumination
+	        do i = 1, n
+		    out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+	    }
+	}
+end
+
+
+# COR2 -- Correct lines for readout axis 2 (columns).  This procedure is
+# more complex than when the readout is along the image lines because the
+# zero level and/or flat field corrections may be single readout column
+# vectors.
+
+procedure cor2s (line, cors, out, overscan, zero, dark, flat, illum,
+	fringe, n, zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+int	line		# Line to be corrected
+int	cors[ARB]	# Correction flags
+short	out[n]		# Output data
+real	overscan[n]	# Overscan value
+short	zero[n]		# Zero level correction
+short	dark[n]		# Dark count correction
+short	flat[n]		# Flat field correction
+short	illum[n]	# Illumination correction
+short	fringe[n]	# Fringe correction
+int	n		# Number of pixels
+pointer	zeroim		# Zero level IMIO pointer (NULL if 1D vector)
+pointer	flatim		# Flat field IMIO pointer (NULL if 1D vector)
+real	darkscale	# Dark count scale factor
+real	flatscale	# Flat field scale factor
+real	illumscale	# Illumination scale factor
+real	frgscale	# Fringe scale factor
+
+short	zeroval
+real	flatval
+int	i, op
+
+begin
+	op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+	switch (op) {
+	case O: # overscan
+	    do i = 1, n
+	        out[i] = out[i] - overscan[i]
+	case Z: # zero level
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - zero[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - zeroval
+	    }
+
+	case ZO: # zero level + overscan
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zero[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zeroval
+	    }
+
+	case D: # dark count
+	    do i = 1, n
+		out[i] = out[i] - darkscale * dark[i]
+	case DO: # dark count + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan[i] - darkscale * dark[i]
+	case DZ: # dark count + zero level
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - zero[i] - darkscale * dark[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - zeroval - darkscale * dark[i]
+	    }
+	case DZO: # dark count + zero level + overscan
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zero[i] -
+			darkscale * dark[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zeroval -
+			darkscale * dark[i]
+	    }
+
+	case F: # flat field
+	    if (flatim != NULL) {
+	        do i = 1, n
+		    out[i] = out[i] * flatscale / flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+		    out[i] = out[i] * flatval
+	    }
+	case FO: # flat field + overscan
+	    if (flatim != NULL) {
+	        do i = 1, n
+	            out[i] = (out[i] - overscan[i]) * flatscale / flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+	            out[i] = (out[i] - overscan[i]) * flatval
+	    }
+	case FZ: # flat field + zero level
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval) * flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i]) * flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval) * flatval
+		}
+	    }
+	case FZO: # flat field + zero level + overscan
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i]) *
+			    flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval) *
+			    flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i]) * flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval) * flatval
+		}
+	    }
+	case FD: # flat field + dark count
+	    if (flatim != NULL) {
+	        do i = 1, n
+		    out[i] = (out[i] - darkscale * dark[i]) * flatscale/flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+		    out[i] = (out[i] - darkscale * dark[i]) * flatval
+	    }
+	case FDO: # flat field + dark count + overscan
+	    if (flatim != NULL) {
+	        do i = 1, n
+		    out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+		        flatscale / flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+		    out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+		        flatval
+	    }
+	case FDZ: # flat field + dark count + zero level
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+			    flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+			    flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+			    flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+			    flatval
+		}
+	    }
+	case FDZO: # flat field + dark count + zero level + overscan
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i] -
+			    darkscale * dark[i]) * flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval -
+			    darkscale * dark[i]) * flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i] -
+			    darkscale * dark[i]) * flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval -
+			    darkscale * dark[i]) * flatval
+		}
+	    }
+	}
+
+	# Often these operations will not be performed so test for no
+	# correction rather than go through the switch.
+
+	op = cors[ILLUMCOR] + cors[FRINGECOR]
+	if (op != 0) {
+	    switch (op) {
+	    case I: # illumination
+	        do i = 1, n
+		    out[i] = out[i] * illumscale / illum[i]
+	    case Q: # fringe
+	        do i = 1, n
+		    out[i] = out[i] - frgscale * fringe[i]
+	    case QI: # fringe + illumination
+	        do i = 1, n
+		    out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+	    }
+	}
+end
+
+# COR1 -- Correct image lines with readout axis 1 (lines).
+
+procedure cor1r (cors, out, overscan, zero, dark, flat, illum,
+	fringe, n, darkscale, flatscale, illumscale, frgscale)
+
+int	cors[ARB]	# Correction flags
+real	out[n]		# Output data
+real	overscan	# Overscan value
+real	zero[n]		# Zero level correction
+real	dark[n]		# Dark count correction
+real	flat[n]		# Flat field correction
+real	illum[n]	# Illumination correction
+real	fringe[n]	# Fringe correction
+int	n		# Number of pixels
+real	darkscale	# Dark count scale factor
+real	flatscale	# Flat field scale factor
+real	illumscale	# Illumination scale factor
+real	frgscale	# Fringe scale factor
+
+int	i, op
+
+begin
+	op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+	switch (op) {
+	case O: # overscan
+	    do i = 1, n
+	        out[i] = out[i] - overscan
+	case Z: # zero level
+	    do i = 1, n
+		out[i] = out[i] - zero[i]
+
+	case ZO: # zero level + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan - zero[i]
+
+	case D: # dark count
+	    do i = 1, n
+		out[i] = out[i] - darkscale * dark[i]
+	case DO: # dark count + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan - darkscale * dark[i]
+	case DZ: # dark count + zero level
+	    do i = 1, n
+		out[i] = out[i] - zero[i] - darkscale * dark[i]
+	case DZO: # dark count + zero level + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan - zero[i] - darkscale * dark[i]
+
+	case F: # flat field
+	    do i = 1, n
+		out[i] = out[i] * flatscale / flat[i]
+	case FO: # flat field + overscan
+	    do i = 1, n
+	        out[i] = (out[i] - overscan) * flatscale / flat[i]
+	case FZ: # flat field + zero level
+	    do i = 1, n
+		out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+	case FZO: # flat field + zero level + overscan
+	    do i = 1, n
+		out[i] = (out[i] - overscan - zero[i]) * flatscale /
+		    flat[i]
+	case FD: # flat field + dark count
+	    do i = 1, n
+		out[i] = (out[i] - darkscale * dark[i]) * flatscale / flat[i]
+	case FDO: # flat field + dark count + overscan
+	    do i = 1, n
+		out[i] = (out[i] - overscan - darkscale * dark[i]) *
+		    flatscale / flat[i]
+	case FDZ: # flat field + dark count + zero level
+	    do i = 1, n
+		out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+		    flatscale / flat[i]
+	case FDZO: # flat field + dark count + zero level + overscan
+	    do i = 1, n
+		out[i] = (out[i] - overscan - zero[i] -
+		    darkscale * dark[i]) * flatscale / flat[i]
+	}
+
+	# Often these operations will not be performed so test for no
+	# correction rather than go through the switch.
+
+	op = cors[ILLUMCOR] + cors[FRINGECOR]
+	if (op != 0) {
+	    switch (op) {
+	    case I: # illumination
+	        do i = 1, n
+		    out[i] = out[i] * illumscale / illum[i]
+	    case Q: # fringe
+	        do i = 1, n
+		    out[i] = out[i] - frgscale * fringe[i]
+	    case QI: # fringe + illumination
+	        do i = 1, n
+		    out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+	    }
+	}
+end
+
+
+# COR2 -- Correct lines for readout axis 2 (columns).  This procedure is
+# more complex than when the readout is along the image lines because the
+# zero level and/or flat field corrections may be single readout column
+# vectors.
+
+procedure cor2r (line, cors, out, overscan, zero, dark, flat, illum,
+	fringe, n, zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+int	line		# Line to be corrected
+int	cors[ARB]	# Correction flags
+real	out[n]		# Output data
+real	overscan[n]	# Overscan value
+real	zero[n]		# Zero level correction
+real	dark[n]		# Dark count correction
+real	flat[n]		# Flat field correction
+real	illum[n]	# Illumination correction
+real	fringe[n]	# Fringe correction
+int	n		# Number of pixels
+pointer	zeroim		# Zero level IMIO pointer (NULL if 1D vector)
+pointer	flatim		# Flat field IMIO pointer (NULL if 1D vector)
+real	darkscale	# Dark count scale factor
+real	flatscale	# Flat field scale factor
+real	illumscale	# Illumination scale factor
+real	frgscale	# Fringe scale factor
+
+real	zeroval
+real	flatval
+int	i, op
+
+begin
+	op = cors[OVERSCAN] + cors[ZEROCOR] + cors[DARKCOR] + cors[FLATCOR]
+	switch (op) {
+	case O: # overscan
+	    do i = 1, n
+	        out[i] = out[i] - overscan[i]
+	case Z: # zero level
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - zero[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - zeroval
+	    }
+
+	case ZO: # zero level + overscan
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zero[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zeroval
+	    }
+
+	case D: # dark count
+	    do i = 1, n
+		out[i] = out[i] - darkscale * dark[i]
+	case DO: # dark count + overscan
+	    do i = 1, n
+		out[i] = out[i] - overscan[i] - darkscale * dark[i]
+	case DZ: # dark count + zero level
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - zero[i] - darkscale * dark[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - zeroval - darkscale * dark[i]
+	    }
+	case DZO: # dark count + zero level + overscan
+	    if (zeroim != NULL)
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zero[i] -
+			darkscale * dark[i]
+	    else {
+		zeroval = zero[line]
+	        do i = 1, n
+		    out[i] = out[i] - overscan[i] - zeroval -
+			darkscale * dark[i]
+	    }
+
+	case F: # flat field
+	    if (flatim != NULL) {
+	        do i = 1, n
+		    out[i] = out[i] * flatscale / flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+		    out[i] = out[i] * flatval
+	    }
+	case FO: # flat field + overscan
+	    if (flatim != NULL) {
+	        do i = 1, n
+	            out[i] = (out[i] - overscan[i]) * flatscale / flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+	            out[i] = (out[i] - overscan[i]) * flatval
+	    }
+	case FZ: # flat field + zero level
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i]) * flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval) * flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i]) * flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval) * flatval
+		}
+	    }
+	case FZO: # flat field + zero level + overscan
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i]) *
+			    flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval) *
+			    flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i]) * flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval) * flatval
+		}
+	    }
+	case FD: # flat field + dark count
+	    if (flatim != NULL) {
+	        do i = 1, n
+		    out[i] = (out[i] - darkscale * dark[i]) * flatscale/flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+		    out[i] = (out[i] - darkscale * dark[i]) * flatval
+	    }
+	case FDO: # flat field + dark count + overscan
+	    if (flatim != NULL) {
+	        do i = 1, n
+		    out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+		        flatscale / flat[i]
+	    } else {
+		flatval = flatscale / flat[line]
+	        do i = 1, n
+		    out[i] = (out[i] - overscan[i] - darkscale * dark[i]) *
+		        flatval
+	    }
+	case FDZ: # flat field + dark count + zero level
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+			    flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+			    flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - zero[i] - darkscale * dark[i]) *
+			    flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - zeroval - darkscale * dark[i]) *
+			    flatval
+		}
+	    }
+	case FDZO: # flat field + dark count + zero level + overscan
+	    if (flatim != NULL) {
+	        if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i] -
+			    darkscale * dark[i]) * flatscale / flat[i]
+		} else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval -
+			    darkscale * dark[i]) * flatscale / flat[i]
+		}
+	    } else {
+		flatval = flatscale / flat[line]
+		if (zeroim != NULL) {
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zero[i] -
+			    darkscale * dark[i]) * flatval
+	        } else {
+		    zeroval = zero[line]
+	            do i = 1, n
+		        out[i] = (out[i] - overscan[i] - zeroval -
+			    darkscale * dark[i]) * flatval
+		}
+	    }
+	}
+
+	# Often these operations will not be performed so test for no
+	# correction rather than go through the switch.
+
+	op = cors[ILLUMCOR] + cors[FRINGECOR]
+	if (op != 0) {
+	    switch (op) {
+	    case I: # illumination
+	        do i = 1, n
+		    out[i] = out[i] * illumscale / illum[i]
+	    case Q: # fringe
+	        do i = 1, n
+		    out[i] = out[i] - frgscale * fringe[i]
+	    case QI: # fringe + illumination
+	        do i = 1, n
+		    out[i] = out[i]*illumscale/illum[i] - frgscale*fringe[i]
+	    }
+	}
+end
+
diff --git a/noao/imred/quadred/src/ccdproc/generic/corinput.x b/noao/imred/quadred/src/ccdproc/generic/corinput.x
new file mode 100644
index 00000000..07afaa41
--- /dev/null
+++ b/noao/imred/quadred/src/ccdproc/generic/corinput.x
@@ -0,0 +1,436 @@
+include	<imhdr.h>
+include	"ccdred.h"
+
+
+# CORINPUT -- Get an input image line, fix the bad pixels, and trim.
+# Return the corrected input line in the output array.
+
+procedure corinputs (in, line, ccd, output, ncols)
+
+pointer	in		# Input IMIO pointer
+int	line		# Corrected output line
+pointer	ccd		# CCD pointer
+short	output[ncols]	# Output data (returned)
+int	ncols		# Number of output columns
+
+int	i, inline
+pointer	inbuf, imgl2s()
+
+begin
+	# Determine the input line in terms of the trimmed output line.
+	if (IN_SEC(ccd) == NULL)
+	    inline = IN_L1(ccd) + line - 1
+	else {
+	    do i = 1, IN_NSEC(ccd) {
+		if (line < OUT_SL1(ccd,i) || line > OUT_SL2(ccd,i))
+		    next
+		inline = IN_SL1(ccd,i) + line - OUT_SL1(ccd,i)
+		break
+	    }
+	}
+
+	# If there are bad lines call a procedure to fix them.  Otherwise
+	# read the image line directly.
+
+	if (NBADLINES(ccd) != 0)
+	    call lfixs (in, inline, Mems[BADLINES(ccd)], IM_LEN(in,1),
+		IM_LEN(in,2), NBADLINES(ccd), inbuf)
+	else
+	    inbuf = imgl2s (in, inline)
+
+	# IF there are bad columns call a procedure to fix them.
+	if (NBADCOLS(ccd) != 0)
+	    call cfixs (inline, Mems[BADCOLS(ccd)], IM_LEN(in,1),
+		IM_LEN(in,2), NBADCOLS(ccd), Mems[inbuf])
+
+	# Move the pixels to the output line.
+	if (IN_SEC(ccd) == NULL)
+	    call amovs (Mems[inbuf+IN_C1(ccd)-OUT_C1(ccd)], output, ncols)
+	else {
+	    do i = 1, IN_NSEC(ccd) {
+		if (inline < IN_SL1(ccd,i) || inline > IN_SL2(ccd,i))
+		    next
+		call amovs (Mems[inbuf+IN_SC1(ccd,i)-OUT_C1(ccd)],
+		    output[OUT_SC1(ccd,i)], OUT_SC2(ccd,i)-OUT_SC1(ccd,i)+1)
+	    }
+	}
+end
+
+
+# CFIX -- Interpolate across bad columns defined in the bad column array.
+
+procedure cfixs (line, badcols, ncols, nlines, nbadcols, data)
+
+int	line				# Line to be fixed
+short	badcols[2, nlines, nbadcols]	# Bad column array
+int	ncols				# Number of columns
+int	nlines				# Number of lines
+int	nbadcols			# Number of bad column regions
+short	data[ncols]			# Data to be fixed
+
+short	val
+real	del
+int	i, j, col1, col2
+
+begin
+	do i = 1, nbadcols {
+	    col1 = badcols[1, line, i]
+	    if (col1 == 0)			# No bad columns
+		return
+	    col2 = badcols[2, line, i]
+	    if (col1 == 1) {			# Bad first column
+		val = data[col2+1]
+		do j = col1, col2
+		    data[j] = val
+	    } else if (col2 == ncols) {		# Bad last column
+		val = data[col1-1]
+		do j = col1, col2
+		    data[j] = val
+	    } else {				# Interpolate
+		del = (data[col2+1] - data[col1-1]) / (col2 - col1 + 2)
+		val = data[col1-1] + del
+		do j = col1, col2
+		    data[j] = val + (j - col1) * del
+	    }
+	}
+end
+
+
+# LFIX -- Get image line and replace bad pixels by interpolation from
+# neighboring lines.  Internal buffers are used to keep the last fixed
+# line and the next good line.  They are allocated with LFIXINIT and
+# freed with LFIXFREE.
+
+procedure lfixs (im, line, badlines, ncols, nlines, nbadlines, data)
+
+pointer	im				# IMIO pointer
+int	line				# Line to be obtained and fixed
+short	badlines[2,nlines,nbadlines]	# Bad line region array
+int	ncols				# Number of columns in image
+int	nlines				# Number of lines in images
+int	nbadlines			# Number of bad line regions
+pointer	data				# Data line pointer (returned)
+
+real	wt1, wt2
+int	i, nextgood, lastgood, col1, col2
+pointer	imgl2s()
+
+pointer	lastbuf, nextbuf
+common	/lfixcom/ lastbuf, nextbuf, lastgood
+
+begin
+	# If this line has bad pixels replace them.  Otherwise just
+	# read the line.
+
+	if (badlines[1, line, 1] != 0) {
+	    # Save the last line which has already been fixed.
+	    if (line != 1)
+		call amovs (Mems[data], Mems[lastbuf], ncols)
+
+	    # Determine the next line with no bad line pixels.  Note that
+	    # this requirement is overly strict since the bad columns
+	    # may not be the same in neighboring lines.
+
+	    nextgood = 0
+	    do i = line+1, nlines {
+		if (badlines[1, i, 1] == 0) {
+		    nextgood = i
+		    break
+		}
+	    }
+
+	    # If the next good line is not the same as previously
+	    # read the data line and store it in a buffer.
+
+	    if ((nextgood != lastgood) && (nextgood != 0)) {
+		data = imgl2s (im, nextgood)
+		call amovs (Mems[data], Mems[nextbuf], ncols)
+		lastgood = nextgood
+	    }
+
+	    # Get the data line.
+	    data = imgl2s (im, line)
+
+	    # Interpolate the bad columns.  At the ends of the image use
+	    # extension otherwise use linear interpolation.
+
+	    if (line == 1) {				# First line is bad
+		do i = 1, nbadlines {
+		    col1 = badlines[1,line,i] - 1
+		    if (col1 == -1)
+			break
+		    col2 = badlines[2,line,i]
+		    call amovs (Mems[nextbuf+col1], Mems[data+col1],
+			col2-col1)
+		}
+	    } else if (nextgood == 0) {			# Last line is bad
+		do i = 1, nbadlines {
+		    col1 = badlines[1,line,i] - 1
+		    if (col1 == -1)
+			break
+		    col2 = badlines[2,line,i]
+		    call amovs (Mems[lastbuf+col1], Mems[data+col1],
+			col2-col1)
+		}
+	    } else {					# Interpolate
+		wt1 = 1. / (nextgood - line + 1)
+		wt2 = 1. - wt1
+		do i = 1, nbadlines {
+		    col1 = badlines[1,line,i] - 1
+		    if (col1 == -1)
+			break
+		    col2 = badlines[2,line,i] - 1
+		    call awsus (Mems[nextbuf+col1], Mems[lastbuf+col1],
+			Mems[data+col1], col2-col1+1, wt1, wt2)
+		}
+	    }
+	} else
+	    data = imgl2s (im, line)
+end
+
+
+# LFIXINIT -- Allocate internal buffers.
+
+procedure lfixinits (im)
+
+pointer	im		# IMIO pointer
+
+int	lastgood
+pointer	lastbuf, nextbuf
+common	/lfixcom/ lastbuf, nextbuf, lastgood
+
+begin
+	call malloc (lastbuf, IM_LEN(im,1), TY_SHORT)
+	call malloc (nextbuf, IM_LEN(im,1), TY_SHORT)
+	lastgood=0
+end
+
+# LFIXFREE -- Free memory when the last line has been obtained.
+
+procedure lfixfrees ()
+
+int	lastgood
+pointer	lastbuf, nextbuf
+common	/lfixcom/ lastbuf, nextbuf, lastgood
+
+begin
+	call mfree (lastbuf, TY_SHORT)
+	call mfree (nextbuf, TY_SHORT)
+end
+
+# CORINPUT -- Get an input image line, fix the bad pixels, and trim.
+# Return the corrected input line in the output array.
+
+procedure corinputr (in, line, ccd, output, ncols)
+
+pointer	in		# Input IMIO pointer
+int	line		# Corrected output line
+pointer	ccd		# CCD pointer
+real	output[ncols]	# Output data (returned)
+int	ncols		# Number of output columns
+
+int	i, inline
+pointer	inbuf, imgl2r()
+
+begin
+	# Determine the input line in terms of the trimmed output line.
+	if (IN_SEC(ccd) == NULL)
+	    inline = IN_L1(ccd) + line - 1
+	else {
+	    do i = 1, IN_NSEC(ccd) {
+		if (line < OUT_SL1(ccd,i) || line > OUT_SL2(ccd,i))
+		    next
+		inline = IN_SL1(ccd,i) + line - OUT_SL1(ccd,i)
+		break
+	    }
+	}
+
+	# If there are bad lines call a procedure to fix them.  Otherwise
+	# read the image line directly.
+
+	if (NBADLINES(ccd) != 0)
+	    call lfixr (in, inline, Mems[BADLINES(ccd)], IM_LEN(in,1),
+		IM_LEN(in,2), NBADLINES(ccd), inbuf)
+	else
+	    inbuf = imgl2r (in, inline)
+
+	# IF there are bad columns call a procedure to fix them.
+	if (NBADCOLS(ccd) != 0)
+	    call cfixr (inline, Mems[BADCOLS(ccd)], IM_LEN(in,1),
+		IM_LEN(in,2), NBADCOLS(ccd), Memr[inbuf])
+
+	# Move the pixels to the output line.
+	if (IN_SEC(ccd) == NULL)
+	    call amovr (Memr[inbuf+IN_C1(ccd)-OUT_C1(ccd)], output, ncols)
+	else {
+	    do i = 1, IN_NSEC(ccd) {
+		if (inline < IN_SL1(ccd,i) || inline > IN_SL2(ccd,i))
+		    next
+		call amovr (Memr[inbuf+IN_SC1(ccd,i)-OUT_C1(ccd)],
+		    output[OUT_SC1(ccd,i)], OUT_SC2(ccd,i)-OUT_SC1(ccd,i)+1)
+	    }
+	}
+end
+
+
+# CFIX -- Interpolate across bad columns defined in the bad column array.
+
+procedure cfixr (line, badcols, ncols, nlines, nbadcols, data)
+
+int	line				# Line to be fixed
+short	badcols[2, nlines, nbadcols]	# Bad column array
+int	ncols				# Number of columns
+int	nlines				# Number of lines
+int	nbadcols			# Number of bad column regions
+real	data[ncols]			# Data to be fixed
+
+real	val
+real	del
+int	i, j, col1, col2
+
+begin
+	do i = 1, nbadcols {
+	    col1 = badcols[1, line, i]
+	    if (col1 == 0)			# No bad columns
+		return
+	    col2 = badcols[2, line, i]
+	    if (col1 == 1) {			# Bad first column
+		val = data[col2+1]
+		do j = col1, col2
+		    data[j] = val
+	    } else if (col2 == ncols) {		# Bad last column
+		val = data[col1-1]
+		do j = col1, col2
+		    data[j] = val
+	    } else {				# Interpolate
+		del = (data[col2+1] - data[col1-1]) / (col2 - col1 + 2)
+		val = data[col1-1] + del
+		do j = col1, col2
+		    data[j] = val + (j - col1) * del
+	    }
+	}
+end
+
+
+# LFIX -- Get image line and replace bad pixels by interpolation from
+# neighboring lines.  Internal buffers are used to keep the last fixed
+# line and the next good line.  They are allocated with LFIXINIT and
+# freed with LFIXFREE.
+
+procedure lfixr (im, line, badlines, ncols, nlines, nbadlines, data)
+
+pointer	im				# IMIO pointer
+int	line				# Line to be obtained and fixed
+short	badlines[2,nlines,nbadlines]	# Bad line region array
+int	ncols				# Number of columns in image
+int	nlines				# Number of lines in images
+int	nbadlines			# Number of bad line regions
+pointer	data				# Data line pointer (returned)
+
+real	wt1, wt2
+int	i, nextgood, lastgood, col1, col2
+pointer	imgl2r()
+
+pointer	lastbuf, nextbuf
+common	/lfixcom/ lastbuf, nextbuf, lastgood
+
+begin
+	# If this line has bad pixels replace them.  Otherwise just
+	# read the line.
+
+	if (badlines[1, line, 1] != 0) {
+	    # Save the last line which has already been fixed.
+	    if (line != 1)
+		call amovr (Memr[data], Memr[lastbuf], ncols)
+
+	    # Determine the next line with no bad line pixels.  Note that
+	    # this requirement is overly strict since the bad columns
+	    # may not be the same in neighboring lines.
+
+	    nextgood = 0
+	    do i = line+1, nlines {
+		if (badlines[1, i, 1] == 0) {
+		    nextgood = i
+		    break
+		}
+	    }
+
+	    # If the next good line is not the same as previously
+	    # read the data line and store it in a buffer.
+
+	    if ((nextgood != lastgood) && (nextgood != 0)) {
+		data = imgl2r (im, nextgood)
+		call amovr (Memr[data], Memr[nextbuf], ncols)
+		lastgood = nextgood
+	    }
+
+	    # Get the data line.
+	    data = imgl2r (im, line)
+
+	    # Interpolate the bad columns.  At the ends of the image use
+	    # extension otherwise use linear interpolation.
+
+	    if (line == 1) {				# First line is bad
+		do i = 1, nbadlines {
+		    col1 = badlines[1,line,i] - 1
+		    if (col1 == -1)
+			break
+		    col2 = badlines[2,line,i]
+		    call amovr (Memr[nextbuf+col1], Memr[data+col1],
+			col2-col1)
+		}
+	    } else if (nextgood == 0) {			# Last line is bad
+		do i = 1, nbadlines {
+		    col1 = badlines[1,line,i] - 1
+		    if (col1 == -1)
+			break
+		    col2 = badlines[2,line,i]
+		    call amovr (Memr[lastbuf+col1], Memr[data+col1],
+			col2-col1)
+		}
+	    } else {					# Interpolate
+		wt1 = 1. / (nextgood - line + 1)
+		wt2 = 1. - wt1
+		do i = 1, nbadlines {
+		    col1 = badlines[1,line,i] - 1
+		    if (col1 == -1)
+			break
+		    col2 = badlines[2,line,i] - 1
+		    call awsur (Memr[nextbuf+col1], Memr[lastbuf+col1],
+			Memr[data+col1], col2-col1+1, wt1, wt2)
+		}
+	    }
+	} else
+	    data = imgl2r (im, line)
+end
+
+
+# LFIXINIT -- Allocate internal buffers.
+
+procedure lfixinitr (im)
+
+pointer	im		# IMIO pointer
+
+int	lastgood
+pointer	lastbuf, nextbuf
+common	/lfixcom/ lastbuf, nextbuf, lastgood
+
+begin
+	call malloc (lastbuf, IM_LEN(im,1), TY_REAL)
+	call malloc (nextbuf, IM_LEN(im,1), TY_REAL)
+	lastgood=0
+end
+
+# LFIXFREE -- Free memory when the last line has been obtained.
+
+procedure lfixfreer ()
+
+int	lastgood
+pointer	lastbuf, nextbuf
+common	/lfixcom/ lastbuf, nextbuf, lastgood
+
+begin
+	call mfree (lastbuf, TY_REAL)
+	call mfree (nextbuf, TY_REAL)
+end
+
diff --git a/noao/imred/quadred/src/ccdproc/generic/mkpkg b/noao/imred/quadred/src/ccdproc/generic/mkpkg
new file mode 100644
index 00000000..0f12b368
--- /dev/null
+++ b/noao/imred/quadred/src/ccdproc/generic/mkpkg
@@ -0,0 +1,12 @@
+# Make CCDRED Package.
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+libpkg.a:
+	cor.x		ccdred.h
+	corinput.x	ccdred.h <imhdr.h>
+	proc.x		ccdred.h <imhdr.h>
+	;
diff --git a/noao/imred/quadred/src/ccdproc/generic/proc.x b/noao/imred/quadred/src/ccdproc/generic/proc.x
new file mode 100644
index 00000000..0251f4f8
--- /dev/null
+++ b/noao/imred/quadred/src/ccdproc/generic/proc.x
@@ -0,0 +1,678 @@
+include	<imhdr.h>
+include	"ccdred.h"
+
+
+.help proc Feb87 noao.imred.ccdred
+.nf ----------------------------------------------------------------------------
+proc -- Process CCD images
+
+These are the main CCD reduction procedures.  There is one for each
+readout axis (lines or columns) and one for short and real image data.
+They apply corrections for bad pixels, overscan levels, zero levels,
+dark counts, flat field response, illumination response, and fringe
+effects.  The image is also trimmed if it was mapped with an image
+section.  The mean value for the output image is computed when the flat
+field or illumination image is processed to form the scale factor for
+these calibrations in order to avoid reading through these image a
+second time.
+
+The processing information and parameters are specified in the CCD
+structure. The processing operations to be performed are specified by
+the correction array CORS in the ccd structure.  There is one array
+element for each operation with indices defined symbolically by macro
+definitions (see ccdred.h); i.e.  FLATCOR.  The value of the array
+element is an integer bit field in which the bit set is the same as the
+array index; i.e element 3 will have the third bit set for an operation
+with array value 2**(3-1)=4.  If an operation is not to be performed
+the bit is not set and the array element has the numeric value zero.
+Note that the addition of several correction elements gives a unique
+bit field describing a combination of operations.  For efficiency the
+most common combinations are implemented as separate units.
+
+The CCD structure also contains the correction or calibration data
+consisting either pointers to data, IMIO pointers for the calibration
+images, and scale factors.
+
+The processing is performed line-by-line.  The procedure CORINPUT is
+called to get an input line.  This procedure trims and fixes bad pixels by
+interpolation.  The output line and lines from the various calibration
+images are read.  The image vectors as well as the overscan vector and
+the scale factors are passed to the procedure COR (which also
+dereferences the pointer data into simple arrays and variables).  That
+procedure does the actual corrections apart from bad pixel
+corrections.
+
+The final optional step is to add each corrected output line to form a
+mean.  This adds efficiency since the operation is done only if desired
+and the output image data is already in memory so there is no I/O
+penalty.
+
+SEE ALSO
+    ccdred.h, cor, fixpix, setfixpix, setoverscan, settrim,
+    setzero, setdark, setflat, setillum, setfringe
+.endhelp ----------------------------------------------------------------------
+
+
+
+# PROC1 -- Process CCD images with readout axis 1 (lines).
+
+procedure proc1s (ccd)
+
+pointer	ccd		# CCD structure
+
+int	i, line, nlin, ncols, nlines, findmean, rep
+int	c1, c2, l1, l2
+real	overscan, darkscale, flatscale, illumscale, frgscale, mean
+short	minrep
+pointer	in, out, zeroim, darkim, flatim, illumim, fringeim
+pointer	outbuf, overscan_vec, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real	asums()
+pointer	imgl2s(), impl2s(), ccd_gls()
+
+begin
+	# Initialize.  If the correction image is 1D then just get the
+	# data once.
+
+	in = IN_IM(ccd)
+	out = OUT_IM(ccd)
+	nlin = IM_LEN(in,2)
+	ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+	nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixinits (in)
+
+	findmean = CORS(ccd, FINDMEAN)
+	if (findmean == YES)
+	    mean = 0.
+	rep = CORS(ccd, MINREP)
+	if (rep == YES)
+	    minrep = MINREPLACE(ccd)
+
+	overscan_vec = OVERSCAN_VEC(ccd)
+
+	if (CORS(ccd, ZEROCOR) == 0) {
+	    zeroim = NULL
+	    zerobuf = 1
+	} else if (IM_LEN(ZERO_IM(ccd),2) == 1) {
+	    zeroim = NULL
+	    zerobuf = ccd_gls (ZERO_IM(ccd), ZERO_C1(ccd), ZERO_C2(ccd), 1)
+	} else
+	    zeroim = ZERO_IM(ccd)
+
+	if (CORS(ccd, DARKCOR) == 0) {
+	    darkim = NULL
+	    darkbuf = 1
+	} else if (IM_LEN(DARK_IM(ccd),2) == 1) {
+	    darkim = NULL
+	    darkbuf = ccd_gls (DARK_IM(ccd), DARK_C1(ccd), DARK_C2(ccd), 1)
+	    darkscale = FLATSCALE(ccd)
+	} else {
+	    darkim = DARK_IM(ccd)
+	    darkscale = DARKSCALE(ccd)
+	}
+
+	if (CORS(ccd, FLATCOR) == 0) {
+	    flatim = NULL
+	    flatbuf = 1
+	} else if (IM_LEN(FLAT_IM(ccd),2) == 1) {
+	    flatim = NULL
+	    flatbuf = ccd_gls (FLAT_IM(ccd), FLAT_C1(ccd), FLAT_C2(ccd), 1)
+	    flatscale = FLATSCALE(ccd)
+	} else {
+	    flatim = FLAT_IM(ccd)
+	    flatscale = FLATSCALE(ccd)
+	}
+
+	if (CORS(ccd, ILLUMCOR) == 0) {
+	    illumim = NULL
+	    illumbuf = 1
+	} else {
+	    illumim = ILLUM_IM(ccd)
+	    illumscale = ILLUMSCALE(ccd)
+	}
+
+	if (CORS(ccd, FRINGECOR) == 0) {
+	    fringeim = NULL
+	    fringebuf = 1
+	} else {
+	    fringeim = FRINGE_IM(ccd)
+	    frgscale = FRINGESCALE(ccd)
+	}
+
+	# For each line read lines from the input.  Procedure CORINPUT
+	# replaces bad pixels by interpolation and applies a trim to the
+	# input.  Get lines from the output image and from the zero level,
+	# dark count, flat field, illumination, and fringe images.
+	# Call COR1 to do the actual pixel corrections.  Finally, add the
+	# output pixels to a sum for computing the mean.
+	# We must copy data outside of the output data section.
+
+	do line = 2 - OUT_L1(ccd), 0
+	    call amovs (
+		Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	do line = 1, nlines {
+	    outbuf = impl2s (out, OUT_L1(ccd)+line-1)
+	    call corinputs (in, line, ccd, Mems[outbuf], IM_LEN(out,1))
+
+	    outbuf = outbuf + OUT_C1(ccd) - 1
+	    if (overscan_vec != NULL)
+		overscan = Memr[overscan_vec+line-1]
+	    if (zeroim != NULL)
+		zerobuf = ccd_gls (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+		    ZERO_L1(ccd)+line-1)
+	    if (darkim != NULL)
+		darkbuf = ccd_gls (darkim, DARK_C1(ccd), DARK_C2(ccd),
+		    DARK_L1(ccd)+line-1)
+	    if (flatim != NULL)
+		flatbuf = ccd_gls (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+		    FLAT_L1(ccd)+line-1)
+	    if (illumim != NULL)
+		illumbuf = ccd_gls (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+		    ILLUM_L1(ccd)+line-1)
+	    if (fringeim != NULL)
+		fringebuf = ccd_gls (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+		    FRINGE_L1(ccd)+line-1)
+
+	    if (OUT_SEC(ccd) == NULL) {
+		call cor1s (CORS(ccd,1), Mems[outbuf],
+		    overscan, Mems[zerobuf], Mems[darkbuf],
+		    Mems[flatbuf], Mems[illumbuf], Mems[fringebuf], ncols,
+		    darkscale, flatscale, illumscale, frgscale)
+	    } else {
+		do i = 1, IN_NSEC(ccd) {
+		    l1 = OUT_SL1(ccd,i)
+		    l2 = OUT_SL2(ccd,i)
+		    if (line < l1 || line > l2)
+			next
+		    c1 = OUT_SC1(ccd,i) - 1
+		    c2 = OUT_SC2(ccd,i) - 1
+		    ncols = c2 - c1 + 1
+		    if (overscan_vec != NULL)
+			overscan = Memr[overscan_vec+(i-1)*nlin+line-l1]
+		    
+		    call cor1s (CORS(ccd,1), Mems[outbuf+c1],
+			overscan, Mems[zerobuf+c1], Mems[darkbuf+c1],
+			Mems[flatbuf+c1], Mems[illumbuf+c1],
+			Mems[fringebuf+c1], ncols,
+			darkscale, flatscale, illumscale, frgscale)
+		}
+	    }
+
+	    if (rep == YES)
+		call amaxks (Mems[outbuf], minrep, Mems[outbuf], ncols)
+	    if (findmean == YES)
+		mean = mean + asums (Mems[outbuf], ncols)
+	}
+
+	do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+	    call amovs (
+		Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	# Compute the mean from the sum of the output pixels.
+	if (findmean == YES)
+	    MEAN(ccd) = mean / ncols / nlines
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixfrees ()
+end
+
+
+# PROC2 -- Process CCD images with readout axis 2 (columns).
+
+procedure proc2s (ccd)
+
+pointer	ccd		# CCD structure
+
+int	line, ncols, nlines, findmean, rep
+real	darkscale, flatscale, illumscale, frgscale, mean
+short	minrep
+pointer	in, out, zeroim, darkim, flatim, illumim, fringeim
+pointer	outbuf, overscan_vec, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real	asums()
+pointer	imgl2s(), impl2s(), imgs2s(), ccd_gls()
+
+begin
+	# Initialize.  If the correction image is 1D then just get the
+	# data once.
+
+	in = IN_IM(ccd)
+	out = OUT_IM(ccd)
+	ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+	nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixinits (in)
+
+	findmean = CORS(ccd, FINDMEAN)
+	if (findmean == YES)
+	    mean = 0.
+	rep = CORS(ccd, MINREP)
+	if (rep == YES)
+	    minrep = MINREPLACE(ccd)
+
+	overscan_vec = OVERSCAN_VEC(ccd)
+
+	if (CORS(ccd, ZEROCOR) == 0) {
+	    zeroim = NULL
+	    zerobuf = 1
+	} else if (IM_LEN(ZERO_IM(ccd),1) == 1) {
+	    zeroim = NULL
+	    zerobuf = imgs2s (ZERO_IM(ccd), 1, 1, ZERO_L1(ccd), ZERO_L2(ccd))
+	} else
+	    zeroim = ZERO_IM(ccd)
+
+	if (CORS(ccd, DARKCOR) == 0) {
+	    darkim = NULL
+	    darkbuf = 1
+	} else if (IM_LEN(DARK_IM(ccd),1) == 1) {
+	    darkim = NULL
+	    darkbuf = imgs2s (DARK_IM(ccd), 1, 1, DARK_L1(ccd), DARK_L2(ccd))
+	    darkscale = DARKSCALE(ccd)
+	} else {
+	    darkim = DARK_IM(ccd)
+	    darkscale = DARKSCALE(ccd)
+	}
+
+	if (CORS(ccd, FLATCOR) == 0) {
+	    flatim = NULL
+	    flatbuf = 1
+	} else if (IM_LEN(FLAT_IM(ccd),1) == 1) {
+	    flatim = NULL
+	    flatbuf = imgs2s (FLAT_IM(ccd), 1, 1, FLAT_L1(ccd), FLAT_L2(ccd))
+	    flatscale = FLATSCALE(ccd)
+	} else {
+	    flatim = FLAT_IM(ccd)
+	    flatscale = FLATSCALE(ccd)
+	}
+
+	if (CORS(ccd, ILLUMCOR) == 0) {
+	    illumim = NULL
+	    illumbuf = 1
+	} else {
+	    illumim = ILLUM_IM(ccd)
+	    illumscale = ILLUMSCALE(ccd)
+	}
+
+	if (CORS(ccd, FRINGECOR) == 0) {
+	    fringeim = NULL
+	    fringebuf = 1
+	} else {
+	    fringeim = FRINGE_IM(ccd)
+	    frgscale = FRINGESCALE(ccd)
+	}
+
+	# For each line read lines from the input.  Procedure CORINPUT
+	# replaces bad pixels by interpolation and applies a trim to the
+	# input.  Get lines from the output image and from the zero level,
+	# dark count, flat field, illumination, and fringe images.
+	# Call COR2 to do the actual pixel corrections.  Finally, add the
+	# output pixels to a sum for computing the mean.
+	# We must copy data outside of the output data section.
+
+	do line = 2 - OUT_L1(ccd), 0
+	    call amovs (
+		Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	do line = 1, nlines {
+	    outbuf = impl2s (out, OUT_L1(ccd)+line-1)
+	    call corinputs (in, line, ccd, Mems[outbuf], IM_LEN(out,1))
+
+	    outbuf = outbuf + OUT_C1(ccd) - 1
+	    if (zeroim != NULL)
+		zerobuf = ccd_gls (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+		    ZERO_L1(ccd)+line-1)
+	    if (darkim != NULL)
+		darkbuf = ccd_gls (darkim, DARK_C1(ccd), DARK_C2(ccd),
+		    DARK_L1(ccd)+line-1)
+	    if (flatim != NULL)
+		flatbuf = ccd_gls (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+		    FLAT_L1(ccd)+line-1)
+	    if (illumim != NULL)
+		illumbuf = ccd_gls (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+		    ILLUM_L1(ccd)+line-1)
+	    if (fringeim != NULL)
+		fringebuf = ccd_gls (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+		    FRINGE_L1(ccd)+line-1)
+
+	    call cor2s (line, CORS(ccd,1), Mems[outbuf],
+		Memr[overscan_vec], Mems[zerobuf], Mems[darkbuf],
+		Mems[flatbuf], Mems[illumbuf], Mems[fringebuf], ncols,
+		zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+	    if (rep == YES)
+		call amaxks (Mems[outbuf], minrep, Mems[outbuf], ncols)
+	    if (findmean == YES)
+		mean = mean + asums (Mems[outbuf], ncols)
+	}
+
+	do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+	    call amovs (
+		Mems[imgl2s(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Mems[impl2s(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	# Compute the mean from the sum of the output pixels.
+	if (findmean == YES)
+	    MEAN(ccd) = mean / ncols / nlines
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixfrees ()
+end
+
+# PROC1 -- Process CCD images with readout axis 1 (lines).
+
+procedure proc1r (ccd)
+
+pointer	ccd		# CCD structure
+
+int	i, line, nlin, ncols, nlines, findmean, rep
+int	c1, c2, l1, l2
+real	overscan, darkscale, flatscale, illumscale, frgscale, mean
+real	minrep
+pointer	in, out, zeroim, darkim, flatim, illumim, fringeim
+pointer	outbuf, overscan_vec, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real	asumr()
+pointer	imgl2r(), impl2r(), ccd_glr()
+
+begin
+	# Initialize.  If the correction image is 1D then just get the
+	# data once.
+
+	in = IN_IM(ccd)
+	out = OUT_IM(ccd)
+	nlin = IM_LEN(in,2)
+	ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+	nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixinitr (in)
+
+	findmean = CORS(ccd, FINDMEAN)
+	if (findmean == YES)
+	    mean = 0.
+	rep = CORS(ccd, MINREP)
+	if (rep == YES)
+	    minrep = MINREPLACE(ccd)
+
+	overscan_vec = OVERSCAN_VEC(ccd)
+
+	if (CORS(ccd, ZEROCOR) == 0) {
+	    zeroim = NULL
+	    zerobuf = 1
+	} else if (IM_LEN(ZERO_IM(ccd),2) == 1) {
+	    zeroim = NULL
+	    zerobuf = ccd_glr (ZERO_IM(ccd), ZERO_C1(ccd), ZERO_C2(ccd), 1)
+	} else
+	    zeroim = ZERO_IM(ccd)
+
+	if (CORS(ccd, DARKCOR) == 0) {
+	    darkim = NULL
+	    darkbuf = 1
+	} else if (IM_LEN(DARK_IM(ccd),2) == 1) {
+	    darkim = NULL
+	    darkbuf = ccd_glr (DARK_IM(ccd), DARK_C1(ccd), DARK_C2(ccd), 1)
+	    darkscale = FLATSCALE(ccd)
+	} else {
+	    darkim = DARK_IM(ccd)
+	    darkscale = DARKSCALE(ccd)
+	}
+
+	if (CORS(ccd, FLATCOR) == 0) {
+	    flatim = NULL
+	    flatbuf = 1
+	} else if (IM_LEN(FLAT_IM(ccd),2) == 1) {
+	    flatim = NULL
+	    flatbuf = ccd_glr (FLAT_IM(ccd), FLAT_C1(ccd), FLAT_C2(ccd), 1)
+	    flatscale = FLATSCALE(ccd)
+	} else {
+	    flatim = FLAT_IM(ccd)
+	    flatscale = FLATSCALE(ccd)
+	}
+
+	if (CORS(ccd, ILLUMCOR) == 0) {
+	    illumim = NULL
+	    illumbuf = 1
+	} else {
+	    illumim = ILLUM_IM(ccd)
+	    illumscale = ILLUMSCALE(ccd)
+	}
+
+	if (CORS(ccd, FRINGECOR) == 0) {
+	    fringeim = NULL
+	    fringebuf = 1
+	} else {
+	    fringeim = FRINGE_IM(ccd)
+	    frgscale = FRINGESCALE(ccd)
+	}
+
+	# For each line read lines from the input.  Procedure CORINPUT
+	# replaces bad pixels by interpolation and applies a trim to the
+	# input.  Get lines from the output image and from the zero level,
+	# dark count, flat field, illumination, and fringe images.
+	# Call COR1 to do the actual pixel corrections.  Finally, add the
+	# output pixels to a sum for computing the mean.
+	# We must copy data outside of the output data section.
+
+	do line = 2 - OUT_L1(ccd), 0
+	    call amovr (
+		Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	do line = 1, nlines {
+	    outbuf = impl2r (out, OUT_L1(ccd)+line-1)
+	    call corinputr (in, line, ccd, Memr[outbuf], IM_LEN(out,1))
+
+	    outbuf = outbuf + OUT_C1(ccd) - 1
+	    if (overscan_vec != NULL)
+		overscan = Memr[overscan_vec+line-1]
+	    if (zeroim != NULL)
+		zerobuf = ccd_glr (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+		    ZERO_L1(ccd)+line-1)
+	    if (darkim != NULL)
+		darkbuf = ccd_glr (darkim, DARK_C1(ccd), DARK_C2(ccd),
+		    DARK_L1(ccd)+line-1)
+	    if (flatim != NULL)
+		flatbuf = ccd_glr (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+		    FLAT_L1(ccd)+line-1)
+	    if (illumim != NULL)
+		illumbuf = ccd_glr (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+		    ILLUM_L1(ccd)+line-1)
+	    if (fringeim != NULL)
+		fringebuf = ccd_glr (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+		    FRINGE_L1(ccd)+line-1)
+
+	    if (OUT_SEC(ccd) == NULL) {
+		call cor1r (CORS(ccd,1), Memr[outbuf],
+		    overscan, Memr[zerobuf], Memr[darkbuf],
+		    Memr[flatbuf], Memr[illumbuf], Memr[fringebuf], ncols,
+		    darkscale, flatscale, illumscale, frgscale)
+	    } else {
+		do i = 1, IN_NSEC(ccd) {
+		    l1 = OUT_SL1(ccd,i)
+		    l2 = OUT_SL2(ccd,i)
+		    if (line < l1 || line > l2)
+			next
+		    c1 = OUT_SC1(ccd,i) - 1
+		    c2 = OUT_SC2(ccd,i) - 1
+		    ncols = c2 - c1 + 1
+		    if (overscan_vec != NULL)
+			overscan = Memr[overscan_vec+(i-1)*nlin+line-l1]
+		    
+		    call cor1r (CORS(ccd,1), Memr[outbuf+c1],
+			overscan, Memr[zerobuf+c1], Memr[darkbuf+c1],
+			Memr[flatbuf+c1], Memr[illumbuf+c1],
+			Memr[fringebuf+c1], ncols,
+			darkscale, flatscale, illumscale, frgscale)
+		}
+	    }
+
+	    if (rep == YES)
+		call amaxkr (Memr[outbuf], minrep, Memr[outbuf], ncols)
+	    if (findmean == YES)
+		mean = mean + asumr (Memr[outbuf], ncols)
+	}
+
+	do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+	    call amovr (
+		Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	# Compute the mean from the sum of the output pixels.
+	if (findmean == YES)
+	    MEAN(ccd) = mean / ncols / nlines
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixfreer ()
+end
+
+
+# PROC2 -- Process CCD images with readout axis 2 (columns).
+
+procedure proc2r (ccd)
+
+pointer	ccd		# CCD structure
+
+int	line, ncols, nlines, findmean, rep
+real	darkscale, flatscale, illumscale, frgscale, mean
+real	minrep
+pointer	in, out, zeroim, darkim, flatim, illumim, fringeim
+pointer	outbuf, overscan_vec, zerobuf, darkbuf, flatbuf, illumbuf, fringebuf
+
+real	asumr()
+pointer	imgl2r(), impl2r(), imgs2r(), ccd_glr()
+
+begin
+	# Initialize.  If the correction image is 1D then just get the
+	# data once.
+
+	in = IN_IM(ccd)
+	out = OUT_IM(ccd)
+	ncols = OUT_C2(ccd) - OUT_C1(ccd) + 1
+	nlines = OUT_L2(ccd) - OUT_L1(ccd) + 1
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixinitr (in)
+
+	findmean = CORS(ccd, FINDMEAN)
+	if (findmean == YES)
+	    mean = 0.
+	rep = CORS(ccd, MINREP)
+	if (rep == YES)
+	    minrep = MINREPLACE(ccd)
+
+	overscan_vec = OVERSCAN_VEC(ccd)
+
+	if (CORS(ccd, ZEROCOR) == 0) {
+	    zeroim = NULL
+	    zerobuf = 1
+	} else if (IM_LEN(ZERO_IM(ccd),1) == 1) {
+	    zeroim = NULL
+	    zerobuf = imgs2r (ZERO_IM(ccd), 1, 1, ZERO_L1(ccd), ZERO_L2(ccd))
+	} else
+	    zeroim = ZERO_IM(ccd)
+
+	if (CORS(ccd, DARKCOR) == 0) {
+	    darkim = NULL
+	    darkbuf = 1
+	} else if (IM_LEN(DARK_IM(ccd),1) == 1) {
+	    darkim = NULL
+	    darkbuf = imgs2r (DARK_IM(ccd), 1, 1, DARK_L1(ccd), DARK_L2(ccd))
+	    darkscale = DARKSCALE(ccd)
+	} else {
+	    darkim = DARK_IM(ccd)
+	    darkscale = DARKSCALE(ccd)
+	}
+
+	if (CORS(ccd, FLATCOR) == 0) {
+	    flatim = NULL
+	    flatbuf = 1
+	} else if (IM_LEN(FLAT_IM(ccd),1) == 1) {
+	    flatim = NULL
+	    flatbuf = imgs2r (FLAT_IM(ccd), 1, 1, FLAT_L1(ccd), FLAT_L2(ccd))
+	    flatscale = FLATSCALE(ccd)
+	} else {
+	    flatim = FLAT_IM(ccd)
+	    flatscale = FLATSCALE(ccd)
+	}
+
+	if (CORS(ccd, ILLUMCOR) == 0) {
+	    illumim = NULL
+	    illumbuf = 1
+	} else {
+	    illumim = ILLUM_IM(ccd)
+	    illumscale = ILLUMSCALE(ccd)
+	}
+
+	if (CORS(ccd, FRINGECOR) == 0) {
+	    fringeim = NULL
+	    fringebuf = 1
+	} else {
+	    fringeim = FRINGE_IM(ccd)
+	    frgscale = FRINGESCALE(ccd)
+	}
+
+	# For each line read lines from the input.  Procedure CORINPUT
+	# replaces bad pixels by interpolation and applies a trim to the
+	# input.  Get lines from the output image and from the zero level,
+	# dark count, flat field, illumination, and fringe images.
+	# Call COR2 to do the actual pixel corrections.  Finally, add the
+	# output pixels to a sum for computing the mean.
+	# We must copy data outside of the output data section.
+
+	do line = 2 - OUT_L1(ccd), 0
+	    call amovr (
+		Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	do line = 1, nlines {
+	    outbuf = impl2r (out, OUT_L1(ccd)+line-1)
+	    call corinputr (in, line, ccd, Memr[outbuf], IM_LEN(out,1))
+
+	    outbuf = outbuf + OUT_C1(ccd) - 1
+	    if (zeroim != NULL)
+		zerobuf = ccd_glr (zeroim, ZERO_C1(ccd), ZERO_C2(ccd),
+		    ZERO_L1(ccd)+line-1)
+	    if (darkim != NULL)
+		darkbuf = ccd_glr (darkim, DARK_C1(ccd), DARK_C2(ccd),
+		    DARK_L1(ccd)+line-1)
+	    if (flatim != NULL)
+		flatbuf = ccd_glr (flatim, FLAT_C1(ccd), FLAT_C2(ccd),
+		    FLAT_L1(ccd)+line-1)
+	    if (illumim != NULL)
+		illumbuf = ccd_glr (illumim, ILLUM_C1(ccd), ILLUM_C2(ccd),
+		    ILLUM_L1(ccd)+line-1)
+	    if (fringeim != NULL)
+		fringebuf = ccd_glr (fringeim, FRINGE_C1(ccd), FRINGE_C2(ccd),
+		    FRINGE_L1(ccd)+line-1)
+
+	    call cor2r (line, CORS(ccd,1), Memr[outbuf],
+		Memr[overscan_vec], Memr[zerobuf], Memr[darkbuf],
+		Memr[flatbuf], Memr[illumbuf], Memr[fringebuf], ncols,
+		zeroim, flatim, darkscale, flatscale, illumscale, frgscale)
+
+	    if (rep == YES)
+		call amaxkr (Memr[outbuf], minrep, Memr[outbuf], ncols)
+	    if (findmean == YES)
+		mean = mean + asumr (Memr[outbuf], ncols)
+	}
+
+	do line = nlines+1, IM_LEN(out,2)-OUT_L1(ccd)+1
+	    call amovr (
+		Memr[imgl2r(in,IN_L1(ccd)+line-1)+IN_C1(ccd)-OUT_C1(ccd)],
+		Memr[impl2r(out,OUT_L1(ccd)+line-1)], IM_LEN(out,1))
+
+	# Compute the mean from the sum of the output pixels.
+	if (findmean == YES)
+	    MEAN(ccd) = mean / ncols / nlines
+
+	if (CORS(ccd, FIXPIX) == YES)
+	    call lfixfreer ()
+end
+