Repatch (from linux) of OSX IRAF

19 files changed, 3046 insertions, 0 deletions

diff --git a/pkg/plot/crtpict/calchgms.x b/pkg/plot/crtpict/calchgms.x
new file mode 100644
index 00000000..15c5aa5a
--- /dev/null
+++ b/pkg/plot/crtpict/calchgms.x
@@ -0,0 +1,192 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include "wdes.h"
+include	"crtpict.h"
+
+# CRT_LINEAR_HGRAM -- Calculate two histograms of an image.  One histogram
+# shows the distribution of intensities in the untransformed image; the other
+# shows the distribution of greyscale values in the transformed image. This
+# procedure assumes a linear transformation.
+
+procedure crt_linear_hgram (im, gp, z1, z2, ztrans, inten_hgram,greys_hgram)
+
+pointer	im		 	# Pointer to image
+pointer	gp			# Graphics descriptor
+real	z1, z2			# Range of intensities mapped
+int	ztrans			# Type of transfer function - linear or unitary
+int	inten_hgram[NBINS]	# Output array of intensity hgram values
+int	greys_hgram[NBINS]	# Output array of greyscale hgram values
+
+pointer buf
+int	npix, nsig_bits, zrange, mask, min_val, max_val
+long	v[IM_MAXDIM]
+int	dz1, dz2, high_zi, low_zi
+real	high_z, low_z
+bool	ggetb()
+pointer	imgnlr(), imgnli()
+int	ggeti()
+errchk 	im_minmax, ggeti, imgnli, imgnlr
+
+begin
+	# If z1 and z2 not in graphcap, set to some reasonable numbers for
+	# plots to be generated.
+	if (ggetb (gp, "z1") && ggetb (gp, "z2")) {
+	    dz1 = ggeti (gp, "z1")
+	    dz2 = ggeti (gp, "z2")
+	} else {
+	    dz1 = 0
+	    dz2 = 255
+	}
+
+	# Calculate number of bits of depth in output device
+	zrange = ggeti (gp, "zr")
+	for (nsig_bits = 0; ; nsig_bits = nsig_bits + 1) {
+	    zrange = zrange / 2
+	    if (zrange == 0)
+		break
+	}
+	mask = (2 ** (nsig_bits)) - 1
+
+	call aclri (inten_hgram, NBINS)
+	call aclri (greys_hgram, NBINS)
+	call amovkl (long(1), v, IM_MAXDIM)
+
+	# Read lines into buffer and accumulate histograms.
+	npix = IM_LEN(im,1)
+
+	if (ztrans == W_UNITARY) {
+	    min_val = int (IM_MIN(im))
+	    max_val = int (IM_MAX(im))
+	    while (imgnli (im, buf, v) != EOF) {
+		call ahgmi (Memi[buf], npix, inten_hgram, NBINS, min_val,
+		    max_val)
+		call aandki  (Memi[buf], mask, Memi[buf], npix)
+		call ahgmi  (Memi[buf], npix, greys_hgram, NBINS, dz1, dz2)
+	    }
+	} else if (IM_PIXTYPE(im) == TY_SHORT) {
+	    min_val = int (IM_MIN(im))
+	    max_val = int (IM_MAX(im))
+	    if (z2 > z1) {
+		# Positive contrast
+		high_zi = int (z2)
+		low_zi  = int (z1)
+	    } else {
+		# Negative contrast
+		high_zi = int (z1)
+		low_zi  = int (z2)
+	    }
+	    while (imgnli (im, buf, v) != EOF) {
+		call ahgmi (Memi[buf], npix, inten_hgram, NBINS, min_val, 
+		    max_val)
+		call ahgmi (Memi[buf], npix, greys_hgram, NBINS, low_zi,
+		    high_zi)
+	    }
+	} else {
+	    if (z2 > z1) {
+		# Positive contrast
+		high_z = z2
+		low_z  = z1
+	    } else {
+		# Negative contrast
+		high_z = z1
+		low_z  = z2
+	    }
+	    while (imgnlr (im, buf, v) != EOF) {
+	        call ahgmr (Memr[buf], npix, inten_hgram, NBINS, IM_MIN(im),
+		    IM_MAX(im))
+	        call ahgmr (Memr[buf], npix, greys_hgram, NBINS, low_z, high_z)
+	    }
+	} 
+end
+
+
+# CRT_USER_HGRAM -- Calculate two histograms of an image.  One histogram
+# shows the distribution of intensities in the untransformed image; the other
+# shows the distribution of greyscale values in the transformed image. This
+# procedure does not assume a linear transformation, but rather uses a user
+# specified look up table.
+
+procedure crt_user_hgram (im, gp, z1, z2, lut, inten_hgram, greys_hgram)
+
+pointer	im		 	# Pointer to image
+pointer	gp			# Graphics descriptor
+real	z1, z2			# Range of intensities mapped
+short	lut[ARB]		# Look up table previously calculated
+int	inten_hgram[NBINS]	# Output array of intensity hgram values
+int	greys_hgram[NBINS]	# Output array of greyscale hgram values
+
+pointer buf, ibuf, sp, rlut
+short	min_val, max_val, short_min, short_max, dz1, dz2
+int	npix
+long	v[IM_MAXDIM]
+short	high_zi, low_zi
+real	high_z, low_z
+pointer	imgnlr(), imgnls()
+errchk 	im_minmax, imgnls, imgnlr
+
+begin
+	# Get max and min in look up table
+	call alims (lut, SZ_BUF, dz1, dz2)
+
+	call aclri (inten_hgram, NBINS)
+	call aclri (greys_hgram, NBINS)
+	call amovkl (long(1), v, IM_MAXDIM)
+
+	# Read lines into buffer and accumulate histograms.
+	npix = IM_LEN(im,1)
+
+	if (IM_PIXTYPE(im) == TY_SHORT) {
+	    min_val = short (IM_MIN(im))
+	    max_val = short (IM_MAX(im))
+	    short_min = short (STARTPT)
+	    short_max = short (ENDPT)
+
+	    if (z2 > z1) {
+		# Positive contrast
+		high_zi = short (z2)
+		low_zi  = short (z1)
+	    } else {
+		# Negative contrast
+		high_zi = short (z1)
+		low_zi  = short (z2)
+	    }
+
+	    while (imgnls (im, buf, v) != EOF) {
+		call ahgms (Mems[buf], npix, inten_hgram, NBINS, min_val, 
+		    max_val)
+		call amaps  (Mems[buf], Mems[buf], npix, low_zi, high_zi, 
+		    short_min, short_max)
+		call aluts  (Mems[buf], Mems[buf], npix, lut)
+		call ahgms (Mems[buf], npix, greys_hgram, NBINS, dz1, dz2)
+	    }
+	} else {
+	    if (z2 > z1) {
+		# Positive contrast
+		high_z = z2
+		low_z  = z1
+	    } else {
+		# Negative contrast
+		high_z = z1
+		low_z  = z2
+	    }
+
+	    call smark (sp)
+	    call salloc (ibuf, npix, TY_INT)
+	    call salloc (rlut, SZ_BUF, TY_REAL)
+	    call achtsr (lut, Memr[rlut], SZ_BUF)
+
+	    while (imgnlr (im, buf, v) != EOF) {
+	        call ahgmr (Memr[buf], npix, inten_hgram, NBINS, IM_MIN(im),
+		    IM_MAX(im))
+
+		call amapr (Memr[buf], Memr[buf], npix, z1, z2, STARTPT, ENDPT)
+		call achtri (Memr[buf], Memi[ibuf], npix)
+		call alutr (Memi[ibuf], Memr[buf], npix, Memr[rlut])
+	        call ahgmr (Memr[buf], npix, greys_hgram, NBINS, real (dz1),
+		    real (dz2))
+	    }
+
+	    call sfree (sp)
+	} 
+end
diff --git a/pkg/plot/crtpict/crtpict.h b/pkg/plot/crtpict/crtpict.h
new file mode 100644
index 00000000..490f45d7
--- /dev/null
+++ b/pkg/plot/crtpict/crtpict.h
@@ -0,0 +1,43 @@
+define	THETA_X			0	# Orientation angle of x axis label
+define	THETA_Y			90	# Orientation angle of y axis label
+define	STEP			10	# Tick marks are placed every 10 pixels
+define	LABEL			100	# Labelled ticks are multiples of 100
+define	SZ_LABEL  		5	# Max number of characters in label
+define	TEXT1			0.15
+define	HGRAMS			0.50
+define	TEXT2			0.10
+define	TEXT3			0.10
+define	SPACE			0.03
+define	NBINS			256	# Number of bins in intensity histogram
+define	SAMPLE_SIZE		1000
+define	STARTPT			0.0E0
+define	ENDPT			4095.0E0
+define	SZ_BUF			4096
+
+define	CRT_XS			0.210
+define	CRT_XE			0.810
+define	CRT_YS			0.076
+define	CRT_YE			0.950
+
+define	LEN_CLPAR 		(15 + 40 + 40 + 20)
+define	FILL			Memi[$1]
+define	REPLICATE		Memi[$1+1]
+define	NSAMPLE_LINES		Memi[$1+2]
+define	LUT			Memi[$1+3]
+define	PERIM			Memi[$1+4]
+define	XMAG			Memr[P2R($1+5)]
+define	YMAG			Memr[P2R($1+6)]
+define	CONTRAST		Memr[P2R($1+7)]
+define	Z1			Memr[P2R($1+8)]
+define	Z2			Memr[P2R($1+9)]
+define	GREYSCALE_FRACTION	Memr[P2R($1+10)]
+define	IMAGE_FRACTION		Memr[P2R($1+11)]
+define	GRAPHICS_FRACTION	Memr[P2R($1+12)]
+define	X_BA			Memr[P2R($1+13)]
+define	Y_BA			Memr[P2R($1+14)]
+define	UFILE			Memc[P2C($1+15)]
+define	ZTRANS			Memc[P2C($1+55)]
+define	DEVICE			Memc[P2C($1+95)]
+
+define	NSTEPS 			16
+define	SZ_LABEL		5
diff --git a/pkg/plot/crtpict/crtpict.semi b/pkg/plot/crtpict/crtpict.semi
new file mode 100644
index 00000000..b2100753
--- /dev/null
+++ b/pkg/plot/crtpict/crtpict.semi
@@ -0,0 +1,263 @@
+# Semicode for the CRTPICT replacement.  Input to the procedure
+# is an IRAF image; output is a file of metacode instructions, 
+# essentially x,y,z for each pixel on the dicomed.  The image intensities
+# are scaled to the dynamic range of the dicomed.  The image
+# is also scaled spatially, either expanded or reduced.  
+
+struct dicomed {
+	
+	# These constants refer to the full dicomed plotting area and will
+	# be read from the graphcap entry for "dicomed".
+
+	int	di_xr	4096	# x resolution
+	int	di_yr	4096	# y resolution
+	int	di_xs	1	# starting x
+	int	di_xe	4096	# ending x
+	int	di_ys	1	# starting y
+	int	di_ye	4096	# ending y
+	int	di_zmin	1	# minimum grey scale value
+	int	di_zmax	254	# maximum grey scale value
+
+	# These constants refer to the dicomed area accessed by crtpict.
+	# They will be stored in file crtpict.h.
+
+	int	crtpict_xr	2059	# x resolution
+	int	crtpict_yr	2931	# y resolution
+	int	crtpict_xs	886	# starting x
+	int	crtpict_xe	2944	# ending x
+	int	crtpict_ys	875	# starting y
+	int	crtpict_ye	3805	# ending y
+}
+
+# The cl parameters that are read_only will be stored in this structure:
+
+struct cl_params {
+
+	bool	fill
+	char	ztrans
+	char	device
+	int	nsample_lines
+	real	xmag
+	real	ymag
+	real	contrast
+	real	z1
+	real	z2
+	real	z1out
+	real	z2out
+	real	greyscale_window
+	real	image_window
+	real	graphics_window
+}
+
+t_crtpict (image_name)
+
+begin
+	cl_params = allocate space for read_only cl_parameters
+
+	# First, get the parameters necessary to calculate z1, z2.
+	if (ztrans = "auto")
+	    get contrast, nsample_lines
+	if (ztrans = "min_max") {
+	    get z1, z2
+	    if (z1 == z2) {
+		get contrast
+		if (abs (contrast) != 1)
+		    get nsample_lines
+	    }
+	}
+
+	# Get parameters necessary to compute the spatial transformation.
+	if (fill) = no {
+	    get xmag, ymag 
+	    if (xmag or ymag < 0)
+		convert them to fractional magnification factors
+	}
+
+	# If the output has been redirected, input is read from the named
+	# command file.  If not, each image name in the input template is
+	# assumed to be preceded by the command "plot".
+	
+	if (output has been redirected) {
+	    redir = true
+	    cmd = open command file
+	}
+
+	# Loop over commands until EOF
+	repeat {
+	    if (redir) {
+		command = get next command from cmd file
+		if (command = EOF)
+		    break
+		if (command != plot) {
+		    reset WCS so gscan coordinates are plotted properly
+		    call gscan (command)
+		} else 
+		    image = image to be plotted
+	    } else {
+		image = next name from expanded template
+		if (image = EOF)
+		    break
+	    }
+
+	    im = open image file
+	    gp = open new graphics descriptor
+
+	    if (user specified output name) {
+	        generate unique output name
+	        call gredir (output name)
+	    }
+		
+	    call plot_image (gp, im, cl_params)
+	}
+
+	close image
+	close gp
+end
+
+
+define	NSTEPS 16
+
+procedure plot_image (gp, im, cl_params)
+
+begin
+	wdes = allocate space for window descriptor as in DISPLAY
+	call establish_transform (gp, im, cl_params, wdes)
+
+	if (cl_params.image_fraction > 0.0)
+	    call transform_image (gp, im, wdes)
+
+	if (cl_params.greyscale_fraction > 0.0)
+	    call draw_greyscale	 (gp, cl_params, NSTEPS)
+
+	if (cl_params.graphics_fraction > 0.0) 
+	    call draw_graphics	 (gp, im, cl_params)
+
+	free (wdes)
+
+
+procedure establish_transform (gp, im, cl_params, wdes)
+
+begin
+	# Procedure xy_scale calculates and stores the spatial
+	# transformation fields of structure wdes
+
+	call xy_scale (gp, cl_params, im, wdes)
+
+	# Now for the intensity to greyscale mapping.  Values z1 and z2,
+	# the intensities that map to the lowest and highest greyscale
+	# levels, are also calculated and stored in the wdes structure. 
+
+	w1 = W_WC(wdes, 1)
+	W_ZT(w1) = W_UNITARY
+
+	if (ztrans = "min_max") {
+	    W_ZT(w1) = W_LINEAR
+	    if (cl_param.z1 != cl_param.z2) {
+		# Use values set by user
+		z1 = cl_param.z1
+		z2 = cl_param.z2
+	    } else {
+		# Use image min and max unless the user has set the contrast
+		# parameter to alter the slope of the transfer function.
+		if (cl_params.contrast == 1) {
+		    z1 = im.im_min
+		    z2 = im.im_max
+		} else if (cl_params.contrast == -1) {
+		        z1 = im.im_max
+		        z2 = im.im_min
+		} else 
+		    call zscale (im, z1, z2, contrast, SAMPLE_SIZE, len_stdline)
+	    }
+	}
+
+	if (ztrans = "auto") {
+	    W_ZT(w1) = W_LINEAR
+	    # Calculate optimum z1 and z2 from image mode
+	    call zscale (im, z1, z2, contrast, SAMPLE_SIZE, len_stdline)
+	}
+
+	W_ZS(w1) = z1
+	W_ZE(w1) = z2
+end
+
+
+procedure xy_scale (gp, cl_params, im, wdes)
+
+begin
+	if (fill) {
+	    # Find the number of device pixels per image pixel required to 
+	    # scale the image to fit the device window.
+	    xscale = scaling factor in x dimension
+	    yscale = scaling factor in y dimension
+
+	    # Preserve the aspect ratio
+	    if (image is longer in x than y)
+		yscale = xscale
+	    else
+		xscale = yscale
+	} else {
+	    # Use the magnification factors specified by the user.
+	    xscale = cl_params.xmag
+	    yscale = cl_params.ymag
+	}
+
+	# The (NDC) device coordinates of the image viewport are stored
+	# as world coordinate system 0.
+	w0 = W_WC(wdes, 0)
+	W_XS(w0) = NDC coord of left edge of viewport
+	W_XE(w0) = NDC coord of right edge of viewport
+	W_YS(w0) = NDC coord of lower edge of viewport
+	W_YE(w0) = NDC coord of upper edge of viewport
+	W_XRES(w0) = number of elements in plotting area x dimension
+	W_YRES(w0) = number of elements in plotting area y dimension
+
+	# The pixel coordinates of the window are stored as world
+	# coordinate system 1.
+	w1 = W_WC(wdes, 1)
+	W_XS(w1) = image column plotted at left edge of window
+	W_XE(w1) = image column plotted at right edge of window
+	W_YS(w1) = image row plotted at lower edge of window
+	W_YE(w1) = image row plotted at upper edge of window
+end
+
+
+procedure draw_greyscale (gp, cl_params, NSTEPS)
+
+begin
+	# The (NDC) device coordinates of the greyscale_window:
+	gs_x1 = crtpict_xs / di_xr
+	gs_x2 = crtpict_x2 / di_xr
+	gs_y1 = im_y2 / di_yr
+	gs_y2 = gs_y1 + ((crtpict_yr * cl_params.greyscale_fraction) / di_yr)
+
+	# Set the viewport and window mapping
+	call gsview (gp, gs_x1, gs_x2, gs_y1, gs_y2)
+	call gswind (gp, 1, NSTEPS, 1, 1)
+
+	# Fill and output greyscale array
+	do i = 1, NSTEPS
+	    grey_levels[i] = grey level 
+
+	call gcell (gp, grey_levels, NSTEPS, 1, 1, 1, NSTEPS, 1)
+end
+
+
+define	NVALUES	256
+
+procedure draw_graphics (gp, im, cl_params)
+
+begin
+	# The (NDC) device coordinates of the graphics viewport:
+	gr_x1 = crtpict_xs / di_xr
+	gr_x2 = crtpict_xe / di_xr
+	gr_y1 = crtpict_ys / di_yr
+	gr_y2 = gr_y1 + ((crtpict_yr * cl_params.graphics_fraction) / di_yr)
+
+	# Set the viewport and window coordinates
+	call gsview (gp, gr_x1, gr_x2, gr_y1, gr_y2)
+	call gswind (gp, 1, crtpict_xr, 1, gr_yr)
+
+	call gtext (for id string, nrows, ncols etc.)
+	call generate_histograms (im, NVALUES)
+	call gploto (?) to plot histograms
+end
diff --git a/pkg/plot/crtpict/crtulut.x b/pkg/plot/crtpict/crtulut.x
new file mode 100644
index 00000000..43609c55
--- /dev/null
+++ b/pkg/plot/crtpict/crtulut.x
@@ -0,0 +1,130 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<error.h>
+include	<ctype.h>
+include	"crtpict.h"
+
+# CRT_ULUT -- Generates a look up table from data supplied by user.  The
+# data is read from a two column text file of intensity, greyscale values.
+# The input data are sorted, then mapped to the x range [0-4096].  A 
+# piecewise linear look up table of 4096 values is then constructed from 
+# the (x,y) pairs given.  A pointer to the look up table, as well as the z1 
+# and z2 intensity endpoints, is returned.
+
+procedure crt_ulut (fname, z1, z2, lut)
+
+char	fname[SZ_FNAME]		# Name of file with intensity, greyscale values
+real	z1			# Intensity mapped to minimum gs value
+real	z2			# Intensity mapped to maximum gs value
+pointer	lut			# Look up table - pointer is returned
+
+pointer	sp, x, y
+int	nvalues, i, j, x1, x2, y1
+real	delta_gs, delta_xv, slope
+errchk	crt_rlut, crt_sort, malloc
+
+begin
+	call smark (sp)
+	call salloc (x, SZ_BUF, TY_REAL)
+	call salloc (y, SZ_BUF, TY_REAL)
+
+	# Read intensities and greyscales from the user's input file.  The
+	# intensity range is then mapped into a standard range and the 
+	# values sorted.
+
+	call crt_rlut (fname, Memr[x], Memr[y], nvalues)
+	call alimr (Memr[x], nvalues, z1, z2)
+	call amapr (Memr[x], Memr[x], nvalues, z1, z2, STARTPT, ENDPT)
+	call crt_sort (Memr[x], Memr[y], nvalues)
+
+	# Fill lut in straight line segments - piecewise linear
+	call malloc (lut, SZ_BUF, TY_SHORT)
+	do i = 1, nvalues-1 {
+	    delta_gs = Memr[y+i] - Memr[y+i-1]
+	    delta_xv = Memr[x+i] - Memr[x+i-1]
+	    slope = delta_gs / delta_xv
+	    x1 = int (Memr[x+i-1]) 
+	    x2 = int (Memr[x+i])
+	    y1 = int (Memr[y+i-1])
+	    do j = x1, x2-1
+		Mems[lut+j-1] = y1 + slope * (j-x1)
+	}
+
+	call sfree (sp)
+end
+
+
+# CRT_RLUT -- Read text file of x, y, values.
+
+procedure crt_rlut (utab, x, y, nvalues)
+
+char	utab[SZ_FNAME]		# Name of list file
+real	x[SZ_BUF]		# Array of x values, filled on return
+real	y[SZ_BUF]		# Array of y values, filled on return
+int	nvalues			# Number of values in x, y vectors - returned
+
+int	n, fd
+pointer	sp, lbuf, ip
+real	xval, yval
+int	getline(), open()
+errchk	open, sscan, getline, malloc
+
+begin
+	call smark (sp)
+	call salloc (lbuf, SZ_LINE, TY_CHAR)
+
+	iferr (fd = open (utab, READ_ONLY, TEXT_FILE))
+	    call error (0, "Error opening user table")
+
+	n = 0
+
+	while (getline (fd, Memc[lbuf]) != EOF) {
+	    # Skip comment lines and blank lines.
+	    if (Memc[lbuf] == '#')
+		next
+	    for (ip=lbuf;  IS_WHITE(Memc[ip]);  ip=ip+1)
+		;
+	    if (Memc[ip] == '\n' || Memc[ip] == EOS)
+		next
+
+	    # Decode the points to be plotted.
+	    call sscan (Memc[ip])
+		call gargr (xval)
+		call gargr (yval)
+
+	    n = n + 1
+	    if (n > SZ_BUF) 
+	        call error (0,
+		    "Intensity transformation table cannot exceed 4096 values")
+
+	    x[n] = xval
+	    y[n] = yval
+	}
+
+	nvalues = n
+	call close (fd)
+	call sfree (sp)
+end
+
+
+# CRT_SORT -- Bubble sort of paired arrays.  
+
+procedure crt_sort (xvals, yvals, nvals)
+
+real	xvals[nvals]		# Array of x values
+real	yvals[nvals]		# Array of y values
+int	nvals			# Number of values in each array
+
+int	i, j
+real	temp
+define	swap	{temp=$1;$1=$2;$2=temp}
+
+begin
+	for (i = nvals; i > 1; i = i - 1)
+	    for (j = 1; j < i; j = j + 1) 
+		if (xvals[j] > xvals[j+1]) {
+		    # Out of order; exchange y values
+		    swap (xvals[j], xvals[j+1])
+		    swap (yvals[j], yvals[j+1])
+		}
+end
diff --git a/pkg/plot/crtpict/drawgraph.x b/pkg/plot/crtpict/drawgraph.x
new file mode 100644
index 00000000..5ee94045
--- /dev/null
+++ b/pkg/plot/crtpict/drawgraph.x
@@ -0,0 +1,153 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<gset.h>
+include	<imhdr.h>
+include	<time.h>
+include "wdes.h"
+include	"crtpict.h"
+
+
+# CRT_DRAW_GRAPHICS -- Draw histogram plots and id information at the bottom of
+# the output print.
+
+procedure crt_draw_graphics (gp, im, cl, wdes)
+
+pointer	gp	# Pointer to graphics descriptor
+pointer	im	# Pointer to input image
+pointer	cl	# Pointer to cl parameter structure
+pointer	wdes	# Pointer to window descriptor
+
+pointer	w1, w0
+char	text[SZ_LINE],  system_id[SZ_LINE]
+int	ndev_rows, ndev_cols
+real	ndc_xs, ndc_xe, ndc_ys, ndc_ye 	 	# Graphics (NDC) viewport
+real	h_xs, h_xe, h_ys, h_ye			# Histogram viewport
+real	tx1_xc, tx1_ly, tx2_xs, tx2_ly, tx3_xs, tx3_ly, tx4_xs, tx4_ly
+real	px1, px2, py1, py2, pxcenter, pycenter, yres
+real	vx1, vx2, vy1, vy2, tx_start, xrf, yrf
+
+int	junk
+pointer	sp, buf
+int	envfind(), envputs()
+
+int	strlen(), ggeti()
+errchk	strlen, ggeti, crt_plot_hgrams, ggwind, ggview
+errchk	gtext
+
+begin
+	ndev_rows = ggeti (gp, "yr")
+	ndev_cols = ggeti (gp, "xr")
+	yres = (CRT_YE - CRT_YS) * ndev_rows
+
+	w0 = W_WC (wdes, 0)
+	w1 = W_WC (wdes, 1)
+
+	# The (NDC) device coordinates of the entire graphics viewport:
+	ndc_xs = CRT_XS 
+	ndc_xe = CRT_XE
+	ndc_ys = CRT_YS
+	ndc_ye = CRT_YS + real (yres * GRAPHICS_FRACTION(cl) / ndev_rows) 
+
+	# Working up from the bottom of the print, locations of various
+	# sections of the graphics are are calculated.  String TEXT1 will
+	# be centered in x at tx1_xc and have lower y coordinate tx1_ly:
+
+	tx1_xc = (ndc_xe + ndc_xs) / 2.0
+	tx1_ly = ndc_ys
+
+	# The three histograms occupy the space calculated next.  This
+	# space is broken into individual plots in a separate procedure.
+
+	h_xs = ndc_xs
+	h_xe = ndc_xe
+	h_ys = ndc_ys + ((ndc_ye - ndc_ys) * (TEXT1 + SPACE))
+	h_ye = h_ys + ((ndc_ye - ndc_ys) * HGRAMS)
+
+	# The left starting position of the text strings is calculated to
+	# line up with the leftmost histogram window:
+	tx_start = h_xs + ((h_xe - h_xs) / 6.0) - ((h_xe - h_xs) / 8.0)
+
+	# String TEXT2 has the following starting_x and lower_y coordinates:
+	tx2_xs = tx_start
+	tx2_ly = h_ye + ((ndc_ye - ndc_ys) * SPACE)
+
+	# String TEXT3 has these starting_x and lower_y coordinates:
+	tx3_xs = tx_start
+	tx3_ly = ndc_ys + ((TEXT1 + HGRAMS + TEXT2 + SPACE) * (ndc_ye - ndc_ys))
+
+	# String TEXT4 has these starting_x and lower_y coordinates:
+	tx4_xs = tx_start
+	tx4_ly = ndc_ys + ((TEXT1+HGRAMS+TEXT2+TEXT3+SPACE) * (ndc_ye - ndc_ys))
+
+	# Draw 3 plots describing transformation of image
+	call crt_plot_histograms (gp, cl, im, wdes, h_xs, h_xe, h_ys, h_ye)
+
+	# Set graphics WCS to WCS 0 for text plotting
+	call gseti (gp, G_WCS, 0)
+
+	# Text line 3 has the image filename and title string.
+	call sprintf (text, SZ_LINE, "%s: %s")
+	    call pargstr (W_IMSECT(wdes))
+	    call pargstr (IM_TITLE(im))
+	call gtext (gp, tx3_xs, tx3_ly, text, "s=0.5")
+
+	# Text line 2 contains image and transformation information; it
+	# is necessary to change to WCS_2 to retrieve the information:
+
+	call gseti (gp, G_WCS, 2)
+	call ggwind (gp, px1, px2, py1, py2)
+	call ggview (gp, vx1, vx2, vy1, vy2)
+	call gseti (gp, G_WCS, 0)
+
+	pxcenter = (px1 + px2) / 2.0
+	pycenter = (py1 + py2) / 2.0
+	xrf = ((vx2 * ndev_cols) - (vx1 * ndev_cols)) / (px2 - px1 + 1.0)
+	yrf = ((vy2 * ndev_rows) - (vy1 * ndev_rows)) / (py2 - py1 + 1.0)
+
+	call sprintf (text, SZ_LINE,
+	    "ncols=%d nrows=%d zmin=%g zmax=%g xc=%0.2f yc=%0.2f")  
+	        call pargi (IM_LEN(im,1))
+	        call pargi (IM_LEN(im,2))
+	        call pargr (IM_MIN(im))
+	        call pargr (IM_MAX(im))
+	        call pargr (pxcenter)
+	        call pargr (pycenter)
+	call sprintf (text[strlen(text)+1], SZ_LINE, " x_rep=%.2f y_rep=%.2f")
+	   call pargr (xrf)
+	   call pargr (yrf)
+	call gtext (gp, tx2_xs, tx2_ly, text, "s=0.35")
+
+	# Text line 1 gives the time and date the output was written
+	call sysid (system_id, SZ_LINE)
+	call gtext (gp, tx1_xc, tx1_ly, system_id, "h=c;s=0.45")
+
+	# Also output transformation information to STDOUT
+	call printf ("ncols=%d nrows=%d zmin=%g zmax=%g xc=%.2f yc=%.2f")
+	    call pargi (IM_LEN(im,1))
+	    call pargi (IM_LEN(im,2))
+	    call pargr (IM_MIN(im))
+	    call pargr (IM_MAX(im))
+	    call pargr (pxcenter)
+	    call pargr (pycenter)
+	call printf (" xrf=%.2f yrf=%.2f\n")
+	    call pargr (xrf)
+	    call pargr (yrf)
+
+	call printf ("%s \n")
+	    call pargstr (system_id)
+
+	# The following was added 17Dec85 at the request of the photo lab.
+	# It allows the negative to be identified easily by user name in
+	# addition to the sequence number written by the 11/23 program.
+
+	call smark (sp)
+	call salloc (buf, SZ_LINE, TY_CHAR)
+
+	if (envfind ("userid", Memc[buf], SZ_LINE) <= 0) {
+	    call getuid (Memc[buf], SZ_LINE)
+	    junk = envputs ("userid", Memc[buf])
+	}
+	call gtext (gp, CRT_XE, 0.001, Memc[buf], "h=r;v=b;s=1.2")
+
+	call sfree (sp)
+end
diff --git a/pkg/plot/crtpict/drawgrey.x b/pkg/plot/crtpict/drawgrey.x
new file mode 100644
index 00000000..6f2807ea
--- /dev/null
+++ b/pkg/plot/crtpict/drawgrey.x
@@ -0,0 +1,63 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<mach.h>
+include	<gset.h>
+include	"wdes.h"
+include	"crtpict.h"
+
+# CRT_DRAW_GREYSCALE -- Draw steps representing greyscale increments on output.
+
+procedure crt_draw_greyscale (gp, cl)
+
+pointer	gp
+pointer	cl
+
+short 	grey[NSTEPS]
+char	label[SZ_LABEL]
+int	ndev_rows, i, dummy
+real	ndc_xs, ndc_xe, ndc_ys, ndc_ye, yres, del_y
+real	delta_grey, delta_x, x_start, x, dz1, dz2
+
+bool	ggetb()
+int	ggeti(), itoc()
+real	ggetr()
+errchk	ggetb, ggeti, ggetr, gpcell, ggetr, gtext
+
+begin
+	if (ggetb (gp, "z1") && ggetb (gp, "z2")) {
+	    dz1 = ggetr (gp, "z1")
+	    dz2 = ggetr (gp, "z2")
+	} else {
+	    dz1 = 0.
+	    dz2 = 255.
+	}
+
+	ndev_rows = ggeti (gp, "yr")
+	yres = (CRT_YE - CRT_YS) * ndev_rows
+
+	# The (NDC) device coordinates of the greyscale_window are calculated.
+
+	ndc_xs = CRT_XS 
+	ndc_xe = CRT_XE
+	ndc_ys = CRT_YS + ((GRAPHICS_FRACTION(cl) + IMAGE_FRACTION(cl) + SPACE)*
+	    yres) / ndev_rows
+	ndc_ye = ndc_ys + ((yres * GREYSCALE_FRACTION(cl)) / ndev_rows) 
+	ndc_ye = min (ndc_ye, CRT_YE)
+	del_y = ndc_ye - ndc_ys
+
+	# Calculate and output grey levels
+	call gseti (gp, G_WCS, 0)
+	delta_grey = (dz2 - dz1) / real(NSTEPS - 1)
+	delta_x = (ndc_xe - ndc_xs) / NSTEPS 
+	x_start = ndc_xs + (delta_x / 2.0)
+	do i = 1, NSTEPS {
+	    grey[i] = short (dz1 + (i-1) * delta_grey + 0.5)
+	    dummy = itoc (int(grey[i]), label, SZ_LABEL)
+	    x = x_start + (i - 1) * delta_x
+	    call gtext (gp, x, ndc_ys, label, "h=c;s=0.25;v=t")
+	}
+	
+	call gpcell (gp, grey, NSTEPS, 1, ndc_xs, ndc_ys + (0.05 * del_y), 
+	    ndc_xe, ndc_ye)
+end
diff --git a/pkg/plot/crtpict/mapimage.x b/pkg/plot/crtpict/mapimage.x
new file mode 100644
index 00000000..9a215278
--- /dev/null
+++ b/pkg/plot/crtpict/mapimage.x
@@ -0,0 +1,172 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<mach.h>
+include	<gset.h>
+include	<error.h>
+include	"wdes.h"
+include	"crtpict.h"
+
+# CRT_MAP_IMAGE -- Output a scaled image window to the device viewport.
+# Spatial scaling is handled by the "scaled input" package, SIGL2[SR]; it is
+# possible to scale the image window to a block averaged device viewport.
+# Image intensities are converted to greyscale values and the input NDC
+# coordinates are "tweaked" to make sure they represent integer device pixels.
+# This tweaking also insures an integer replication factor between image
+# pixels and device pixels.  Type short pixels are treated as a special
+# case to minimize vector operations.
+
+procedure crt_map_image (im, gp, px1,px2,py1,py2, ndc_xs,ndc_xe,ndc_ys,ndc_ye, 
+        nx_output, ny_output, z1,z2,zt, cl)
+
+pointer	im				# input image
+pointer	gp				# graphics descriptor
+real	px1,px2,py1,py2			# input section
+real	ndc_xs,ndc_xe,ndc_ys,ndc_ye	# NDC of output section
+int	nx_output, ny_output		# Number of output pixels.  Image pixels
+					# are scaled to these dimensions.
+real	z1,z2				# range of intensities to be mapped.
+int	zt				# specified greyscale transform type
+pointer	cl				# Pointer to crtpict structure
+
+bool	unitary_greyscale_transformation
+pointer	in, si, sline, rline, llut, sp
+short	sz1, sz2, sdz1, sdz2, lut1, lut2
+real	dz1, dz2, y1, y2, delta_y
+int	ndev_cols, ndev_rows, nline, ny_device
+int	xblk, yblk
+bool	ggetb(), fp_equalr()
+int	ggeti()
+real	ggetr()
+pointer	sigl2s(), sigl2r(), sigl2_setup()
+errchk	sigl2s, sigl2r, sigl2_setup, ndc_tweak_ndc, ggeti, malloc, ggetr
+errchk	ggetb, ggetr, gpcell, crt_ulut
+
+begin
+	call smark (sp)
+	call salloc (sline, nx_output, TY_SHORT)
+	if (IM_PIXTYPE(im) != TY_SHORT)
+	    call salloc (rline, nx_output, TY_REAL)
+
+	# Calculate and allocate heap space needed for an image row.
+	ndev_cols = ggeti (gp, "xr")
+	ndev_rows = ggeti (gp, "yr")
+	ny_device  = ((ndc_ye * ndev_rows) - (ndc_ys * ndev_rows)) + 1
+
+	# This sets up for the scaled image input 
+	xblk = INDEFI
+	yblk = INDEFI
+	si = sigl2_setup (im, px1,px2,nx_output,xblk, py1,py2,ny_output,yblk)
+
+	# If user has supplied look up table, it has to be dealt with at
+	# this point.  Greyscale transform is coming up, and the transfer
+	# function will be plotted at a later step.
+
+	if (zt == W_USER) {
+	    iferr (call crt_ulut (UFILE(cl), z1, z2, llut))
+		call erract (EA_FATAL)
+	    LUT(cl) = llut
+	    call alims (Mems[llut], SZ_BUF, lut1, lut2)
+	}
+
+	# If device can't output greyscale information, return at this point.
+	if (! ggetb (gp, "zr")) {
+	    call eprintf ("Graphics device doesn't support greyscale output\n")
+	    return
+	}
+
+	# Determine the device range for the greyscale transformation.
+	if (ggetb (gp, "z1") && ggetb (gp, "z2")) {
+	    dz1 = ggetr (gp, "z1")
+	    dz2 = ggetr (gp, "z2")
+	} else {
+	    dz1 = 0.
+	    dz2 = 255.
+	}
+
+	# And now a quick test to make sure user specified greyscale and
+	# intensity ranges are reasonable.
+	if (zt == W_USER) {
+	    sdz1 = short (dz1)
+	    sdz2 = short (dz2)
+	    if (lut2 < sdz1 || lut1 > sdz2)
+		call eprintf ("User specified greyscales out of range\n")
+	    if (z2 < IM_MIN(im) || z1 > IM_MAX(im))
+		call eprintf ("User specified intensities out of range\n")
+	}
+
+	if (zt == W_UNITARY)
+	    unitary_greyscale_transformation = true
+	else
+	    unitary_greyscale_transformation =
+		(fp_equalr (dz1,z1) && fp_equalr (dz2,z2)) || fp_equalr (z1,z2)
+
+	# Calculate the delta_y, that is, the change in ndc coordinate
+	# with each output row.  It has been assurred by tweak_ndc that
+	# the ratio of device rows to output pixels is an integer.  
+
+	delta_y = (real (ny_device) / ny_output) / real (ndev_rows)
+
+	# For TY_SHORT pixels, pixel intensities are converted to greyscale
+	# values, then output with gpcell.
+	if (IM_PIXTYPE(im) == TY_SHORT) {
+	    for (nline=1;  nline <= ny_output;  nline=nline+1) {
+		in  = sigl2s (si, nline)
+		if (unitary_greyscale_transformation)
+		    call amovs (Mems[in], Mems[sline], nx_output)
+		else if (zt == W_LINEAR) {
+		    sz1  = short (z1)
+		    sz2  = short (z2)
+		    sdz1 = short (dz1)
+		    sdz2 = short (dz2)
+		    call amaps (Mems[in], Mems[sline], nx_output, sz1, sz2, 
+			sdz1, sdz2)
+		} else if (zt == W_USER) {
+		    sz1  = short (z1)
+		    sz2  = short (z2)
+		    sdz1 = short (STARTPT)
+		    sdz2 = short (ENDPT)
+		    call amaps (Mems[in], Mems[sline], nx_output, sz1, sz2,
+			sdz1, sdz2)
+		    call aluts (Mems[sline], Mems[sline], nx_output, Mems[llut])
+		}
+		
+		# Now put line out to greyscale device
+		y1 = ndc_ys + (nline - 1) * delta_y
+		y2 = ndc_ys + (nline * delta_y)
+
+		call gpcell (gp, Mems[sline], nx_output, 1, ndc_xs, y1, ndc_xe,
+		    y2)
+	    }
+	} else {
+	    # Pixels are treated as TY_REAL; intensities are converted to
+	    # greyscale values, then converted to TY_SHORT for gpcell output.
+	    for (nline=1;  nline <= ny_output;  nline=nline+1) {
+		in  = sigl2r (si, nline)
+		if (unitary_greyscale_transformation) {
+		    call amovr (Memr[in], Memr[rline], nx_output)
+		    call achtrs (Memr[rline], Mems[sline], nx_output)
+		} else if (zt == W_LINEAR) {
+		    call amapr (Memr[in], Memr[rline], nx_output, z1, z2, dz1, 
+			dz2)
+		    call achtrs (Memr[rline], Mems[sline], nx_output)
+		} else if (zt == W_USER) {
+		    call amapr (Memr[in], Memr[rline], nx_output, z1, z2, 
+			STARTPT, ENDPT)
+		    call achtrs (Memr[rline], Mems[sline], nx_output)
+		    call aluts (Mems[sline], Mems[sline], nx_output, Mems[llut])
+		}
+		
+		# Output line to greyscale device
+		y1 = ndc_ys + (nline - 1) * delta_y
+		y2 = ndc_ys + (nline * delta_y)
+
+		call gpcell (gp, Mems[sline], nx_output, 1, ndc_xs, y1, 
+		    ndc_xe, y2)
+	    }
+	}
+
+	# Free allocate memory
+	call sigl2_free (si)
+	call sfree (sp)
+end
diff --git a/pkg/plot/crtpict/minmax.x b/pkg/plot/crtpict/minmax.x
new file mode 100644
index 00000000..092d3b9e
--- /dev/null
+++ b/pkg/plot/crtpict/minmax.x
@@ -0,0 +1,75 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# IM_MINMAX -- Compute the minimum and maximum pixel values of an image.
+# Works for images of any dimensionality, size, or datatype, although
+# the min and max values can currently only be stored in the image header
+# as real values.
+
+procedure im_minmax (im, min_value, max_value)
+
+pointer	im				# image descriptor
+real	min_value			# minimum pixel value in image (out)
+real	max_value			# maximum pixel value in image (out)
+
+pointer	buf
+bool	first_line
+long	v[IM_MAXDIM]
+short	minval_s, maxval_s
+long	minval_l, maxval_l
+real	minval_r, maxval_r
+int	imgnls(), imgnll(), imgnlr()
+errchk	amovkl, imgnls, imgnll, imgnlr, alims, aliml, alimr
+
+begin
+	call amovkl (long(1), v, IM_MAXDIM)		# start vector
+	first_line = true
+	min_value = INDEF
+	max_value = INDEF
+
+	switch (IM_PIXTYPE(im)) {
+	case TY_SHORT:
+	    while (imgnls (im, buf, v) != EOF) {
+		call alims (Mems[buf], IM_LEN(im,1), minval_s, maxval_s)
+		if (first_line) {
+		    min_value = minval_s
+		    max_value = maxval_s
+		    first_line = false
+		} else {
+		    if (minval_s < min_value)
+			min_value = minval_s
+		    if (maxval_s > max_value)
+			max_value = maxval_s
+		}
+	    }
+	case TY_USHORT, TY_INT, TY_LONG:
+	    while (imgnll (im, buf, v) != EOF) {
+		call aliml (Meml[buf], IM_LEN(im,1), minval_l, maxval_l)
+		if (first_line) {
+		    min_value = minval_l
+		    max_value = maxval_l
+		    first_line = false
+		} else {
+		    if (minval_l < min_value)
+			min_value = minval_l
+		    if (maxval_l > max_value)
+			max_value = maxval_l
+		}
+	    }
+	default:
+	    while (imgnlr (im, buf, v) != EOF) {
+		call alimr (Memr[buf], IM_LEN(im,1), minval_r, maxval_r)
+		if (first_line) {
+		    min_value = minval_r
+		    max_value = maxval_r
+		    first_line = false
+		} else {
+		    if (minval_r < min_value)
+			min_value = minval_r
+		    if (maxval_r > max_value)
+			max_value = maxval_r
+		}
+	    }
+	}
+end
diff --git a/pkg/plot/crtpict/mkpkg b/pkg/plot/crtpict/mkpkg
new file mode 100644
index 00000000..9f41a011
--- /dev/null
+++ b/pkg/plot/crtpict/mkpkg
@@ -0,0 +1,24 @@
+# Makelib file for CRTPICT contributions to the plot package library.
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+libpkg.a:
+	calchgms.x	crtpict.h <imhdr.h> wdes.h
+	crtulut.x	crtpict.h <ctype.h> <error.h>
+	drawgraph.x	crtpict.h <gset.h> <imhdr.h> <time.h> wdes.h
+	drawgrey.x	crtpict.h wdes.h <gset.h> <imhdr.h> <mach.h>
+	mapimage.x	crtpict.h wdes.h <error.h> <gset.h> <imhdr.h> <mach.h>
+	minmax.x	<imhdr.h>
+	plothgms.x	crtpict.h <gset.h> <imhdr.h> <mach.h> wdes.h
+	plotimage.x	crtpict.h wdes.h <mach.h>
+	setxform.x	crtpict.h wdes.h <imhdr.h> <mach.h>
+	sigl2.x		<error.h> <imhdr.h>
+	t_crtpict.x	crtpict.h <fset.h> <gset.h> <imhdr.h> <mach.h> <error.h>
+	tweakndc.x	
+	xformimage.x	crtpict.h wdes.h <gset.h> <imhdr.h> <mach.h>
+	xyscale.x	crtpict.h wdes.h <error.h> <imhdr.h> <mach.h>
+	zscale.x	<imhdr.h>
+	;
diff --git a/pkg/plot/crtpict/plothgms.x b/pkg/plot/crtpict/plothgms.x
new file mode 100644
index 00000000..803c709f
--- /dev/null
+++ b/pkg/plot/crtpict/plothgms.x
@@ -0,0 +1,209 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<gset.h>
+include	<imhdr.h>
+include	<mach.h>
+include "wdes.h"
+include	"crtpict.h"
+
+# CRT_PLOT_HISTOGRAMS -- Calculate and plot three histograms describing the
+# intensity to greyscale mapping.
+
+procedure crt_plot_histograms (gp, cl, im, wdes, xs, xe, ys, ye)
+
+pointer	gp
+pointer	cl		# Pointer to cl structure.
+pointer im
+pointer	wdes
+real	xs, xe, ys, ye, z1, z2
+
+pointer	w1, inten_hgram, greys_hgram, sp, text, syv, greys, hgram
+pointer real_greys, xval, yval
+int	nsig_bits, i, zrange, mask
+real	plot_ys, plot_ye, plot_width, plot_spacing, x, y, delta_inten, inten
+real	plot1_xs, plot1_xe, plot2_xs, plot2_xe, plot3_xs, plot3_xe 
+real	dz1, dz2, gio_char_x
+real	major_length, minor_length, ux1, ux2, uy1, uy2, y_pos, label_y
+real	wx1, wx2, wy1, wy2
+
+bool	ggetb()
+real 	ggetr()
+int	ggeti(), and()
+errchk	ggetb, ggeti, ggetr, crt_calc_hgrams, ggwind, gswind, gsview, gploto
+errchk	gsetr, gseti, glabax, amapr, gpline, gvline, achtir, gtext
+
+begin
+	call smark (sp)
+	call salloc (inten_hgram, NBINS, TY_INT)
+	call salloc (greys_hgram, NBINS, TY_INT)
+	call salloc (text, SZ_LINE, TY_CHAR)
+	call salloc (syv, NBINS, TY_SHORT)
+	call salloc (greys, NBINS, TY_INT)
+	call salloc (hgram, NBINS, TY_REAL)
+	call salloc (real_greys, NBINS, TY_REAL)
+	call salloc (xval, NBINS, TY_REAL)
+	call salloc (yval, NBINS, TY_REAL)
+
+	# First, get pointer to WCS 1 and some device parameters.
+	w1 = W_WC(wdes, 1)
+	z1 = W_ZS(w1)
+	z2 = W_ZE(w1)
+
+	# If z1 and z2 not in graphcap, set to some reasonable numbers for
+	# plots to be generated.
+
+	if (ggetb (gp, "z1") && ggetb (gp, "z2")) {
+	    dz1 = ggetr (gp, "z1")
+	    dz2 = ggetr (gp, "z2")
+	} else {
+	    dz1 = 0.
+	    dz2 = 255.
+	}
+
+	# To allow room for annotation, the y limits of each plot are
+	# drawn in by 5%.  The y limits are the same for each plot.
+
+	plot_ys = ys + (0.05 * (ye - ys))
+	plot_ye = ye - (0.05 * (ye - ys))
+	label_y = plot_ys - (plot_ye - plot_ys) * 0.20
+
+	# Now calculate the x limits.  Each plot occupys a fourth of the
+	# available space in x.  The distance between plot centers is a
+	# third of the available space.
+
+	plot_width = (xe - xs) / 4.0
+	plot_spacing = (xe - xs) / 3.0
+
+	plot1_xs = xs + (plot_spacing / 2.0) - (plot_width / 2.0)
+	plot1_xe = plot1_xs + plot_width
+	plot2_xs = plot1_xs + plot_spacing
+	plot2_xe = plot2_xs + plot_width
+	plot3_xs = plot2_xs + plot_spacing
+	plot3_xe = plot3_xs + plot_width
+
+	# Calculate the histograms for both the untransformed (intensity) and
+	# transformed (greyscale) image in a single procedure. A separate
+	# path is taken for linear or user transformations:
+
+	if (W_ZT(w1) == W_USER)
+	    call crt_user_hgram (im, gp, z1, z2, Mems[LUT(cl)], 
+		Memi[inten_hgram], Memi[greys_hgram])
+	else
+	    call crt_linear_hgram (im, gp, z1, z2, W_ZT(w1), Memi[inten_hgram], 
+	        Memi[greys_hgram])
+
+	# Each histogram plot is a separate mapping in WCS 3
+	call gseti (gp, G_WCS, 3)
+
+	# The first histogram shows the number of pixels at a given
+	# intensity versus intensity for the original image.
+
+	call gsview (gp, plot1_xs, plot1_xe, plot_ys, plot_ye)
+	gio_char_x = ((plot1_xe - plot1_xs) / 50.) / ggetr (gp, "cw")
+	major_length = gio_char_x * (ggetr (gp, "cw"))
+	minor_length = gio_char_x * (0.5 * (ggetr (gp, "cw")))
+
+	call gseti (gp, G_YTRAN, GW_LOG)
+	call gseti (gp, G_ROUND, YES)
+	call gsetr (gp, G_MAJORLENGTH, major_length)
+	call gsetr (gp, G_MINORLENGTH, minor_length)
+	call gseti (gp, G_LABELAXIS, NO)
+	call gsetr (gp, G_TICKLABELSIZE, 0.25)
+	call achtir (Memi[inten_hgram], Memr[hgram], NBINS)
+	call gploto (gp, Memr[hgram], NBINS, IM_MIN(im), IM_MAX(im), "")
+
+	# Now to label the plot:
+	call ggwind (gp, ux1, ux2, uy1, uy2)
+	y_pos = uy1 - ((uy2 - uy1) * 0.20)  #y_pos below yaxis by 20% of height
+	call gseti (gp, G_YTRAN, GW_LINEAR)
+	call gtext (gp, (ux1+ux2)/2.0, y_pos, "LOG10(N(DN)) VS DN", 
+	    "v=t;h=c;s=.25")
+
+	# The third plot shows the number of pixels at a given greyscale
+	# level versus greyscale level for the range of intensities
+	# transformed. 
+
+	call gsview (gp, plot3_xs, plot3_xe, plot_ys, plot_ye)
+	call achtir (Memi[greys_hgram], Memr[hgram], NBINS)
+
+	if (z2 > z1)
+	    call gploto (gp, Memr[hgram], NBINS, dz1, dz2, "")
+	else
+	    call gploto (gp, Memr[hgram], NBINS, dz2, dz1, "")
+
+	# Now to label the plot:
+	call ggwind (gp, ux1, ux2, uy1, uy2)
+	call gseti (gp, G_YTRAN, GW_LINEAR)
+	y_pos = uy1 - ((uy2 - uy1) * 0.20) # y_pos below yaxis by 20% of height
+	call gtext (gp, (ux1+ux2)/2.0, y_pos, "TRANSFORMED HISTOGRAM", 
+	    "v=t;h=c;s=.25")
+
+	# The second plot shows how the dynamic range of the transformed
+	# image maps to the dynamic range of the output device.
+
+	call gsview (gp, plot2_xs, plot2_xe, plot_ys, plot_ye)
+	call gswind (gp, IM_MIN(im), IM_MAX(im), real (dz1), real (dz2))
+	call gseti (gp, G_YTRAN, GW_LINEAR)
+	call glabax (gp, "", "", "")
+
+	if (W_ZT(w1) != W_UNITARY) {
+	    do i = 1, NBINS 
+		Memr[xval+i-1] = IM_MIN(im) + (i-1) * (IM_MAX(im) - 
+		    IM_MIN(im))/ (NBINS-1)
+
+	    if (W_ZT(w1) == W_USER) {
+	        call sprintf (Memc[text], SZ_LINE, 
+		    "USER DEFINED FUNCTION: FROM %g TO %g")
+		call pargr (z1)
+		call pargr (z2)
+		call amapr (Memr[xval], Memr[yval], NBINS, z1, z2, STARTPT, 
+		    ENDPT)
+		call achtrs (Memr[yval], Mems[syv], NBINS)
+		call aluts (Mems[syv], Mems[syv], NBINS, Mems[LUT[cl]])
+		call achtsr (Mems[syv], Memr[yval], NBINS)
+	    } else {
+	        call sprintf (Memc[text], SZ_LINE, 
+		    "TRANSFER FUNCTION: LINEAR FROM %g TO %g")
+		call pargr (z1)
+		call pargr (z2)
+	        call amapr (Memr[xval], Memr[yval], NBINS, z1, z2, real (dz1), 
+		    real (dz2))
+	    }
+
+	    call gpline (gp, Memr[xval], Memr[yval], NBINS)
+	    call ggwind (gp, wx1, wx2, wy1, wy2)
+	    x = (wx2 + wx1) / 2.0
+	    y = wy1 - (wy2 - wy1) * 0.20
+	    call gtext (gp, x, y, Memc[text], "h=c;v=t;s=0.25")
+	    call printf ("%s\n")
+		call pargstr (Memc[text])
+
+	} else {
+	    # Calculate number of bits depth in output device
+	    zrange = ggeti (gp, "zr")
+	    for (nsig_bits = 0; ; nsig_bits = nsig_bits + 1) {
+	        zrange = zrange / 2
+	        if (zrange == 0)
+	            break
+	    }
+
+	    # Truncate intensities to dynamic range of output device.
+	    delta_inten = (IM_MAX(im) - IM_MIN(im)) / (NBINS - 1)
+	    mask = 2**(nsig_bits) - 1
+	    do i = 1, NBINS {
+		inten =  IM_MIN(im) + ((i-1) * delta_inten)
+		Memi[greys+i-1] = and (int (inten), mask)
+	    }
+
+	    call achtir (Memi[greys], Memr[real_greys], NBINS)
+	    call gvline (gp, Memr[real_greys], NBINS, IM_MIN(im), IM_MAX(im))
+	    call ggwind (gp, wx1, wx2, wy1, wy2)
+	    x = (wx2 + wx1) / 2.0
+	    y = wy1 - (wy2 - wy1) * 0.20
+	    call gtext (gp, x, y, "TRANSFER FUNCTION: UNITARY","h=c;v=t;s=0.25")
+	    call printf ("Unitary Transfer Function; Lowest %d bits output.\n")
+		call pargi (nsig_bits)
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/plot/crtpict/plotimage.x b/pkg/plot/crtpict/plotimage.x
new file mode 100644
index 00000000..add1ed8a
--- /dev/null
+++ b/pkg/plot/crtpict/plotimage.x
@@ -0,0 +1,40 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	"wdes.h"
+include	"crtpict.h"
+
+# CRT_PLOT_IMAGE - Plot the image, graphics and greyscale portion of
+# each image to be transformed.
+
+procedure crt_plot_image (gp, im, image, cl)
+
+pointer	gp			# Graphics descriptor
+pointer	im			# Pointer to image
+char	image[SZ_FNAME]		# Image filename
+pointer	cl			# Pointer to structure of cl parameters
+
+pointer	sp, wdes
+errchk	crt_establish_transform, crt_transform_image
+errchk  crt_draw_graphics, crt_draw_greyscale
+
+begin
+	call smark (sp)
+	call salloc (wdes, LEN_WDES, TY_STRUCT)
+	call strcpy (image, W_IMSECT(wdes), W_SZIMSECT)
+
+	if (IMAGE_FRACTION(cl) > EPSILON) {
+	    call crt_establish_transform (gp, im, cl, wdes)
+	    call crt_transform_image (gp, im, wdes, cl)
+	}
+
+	if (GRAPHICS_FRACTION(cl) > EPSILON)  {
+	    call crt_draw_graphics (gp, im, cl, wdes)
+	}
+
+	if (GREYSCALE_FRACTION(cl) > EPSILON) {
+	    call crt_draw_greyscale (gp, cl)
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/plot/crtpict/setxform.x b/pkg/plot/crtpict/setxform.x
new file mode 100644
index 00000000..a6c50dfb
--- /dev/null
+++ b/pkg/plot/crtpict/setxform.x
@@ -0,0 +1,96 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<mach.h>
+include	"wdes.h"
+include	"crtpict.h"
+
+
+# CRT_ESTABLISH_TRANSFORM -- Set up both the spatial and greyscale mapping.
+# The window descriptor "wdes" is filled.
+
+procedure crt_establish_transform (gp, im, cl, wdes)
+
+pointer	gp
+pointer im
+pointer	cl
+pointer	wdes
+
+int	w1, len_stdline
+real	z1, z2
+
+bool	fp_equalr()
+int	strncmp()
+errchk	crt_xy_scale, strcmp, im_minmax, zscale
+
+begin
+	# Procedure xy_scale calculates and stores the spatial
+	# transformation fields of structure wdes.
+
+	call crt_xy_scale (gp, cl, im, wdes)
+
+	w1 = W_WC(wdes, 1)
+
+	# Now for the intensity to greyscale mapping.  Values z1 and z2,
+	# the intensities that map to the lowest and highest greyscale
+	# levels, are calculated and stored in the wdes structure. First,
+	# set up the default values.  These will not be changed if 
+	# ztrans = "user".
+
+	W_ZT(w1) = W_USER
+	z1 = INDEFR
+	z2 = INDEFR
+
+	# Put up to date min, max values in the header structure, if necessary
+	if (IM_LIMTIME(im) < IM_MTIME(im))
+	    call im_minmax (im, IM_MIN(im), IM_MAX(im))
+
+	if (strncmp (ZTRANS(cl), "min_max", 1) == 0) {
+	    W_ZT(w1) = W_LINEAR
+	    if (Z1(cl) != Z2(cl)) {
+		# Use values set by user
+		z1 = Z1(cl)
+		z2 = Z2(cl)
+	    } else {
+		# Use image min and max unless the user has set the contrast
+		# parameter to alter the slope of the transfer function.
+		if (abs (CONTRAST(cl) - 1.0) > EPSILON) {
+		    # CONTRAST = 1.0
+		    z1 = IM_MIN(im)
+		    z2 = IM_MAX(im)
+		} else if (abs (CONTRAST(cl) + 1.0) > EPSILON) {
+		    # CONTRAST = -1.0
+		    z1 = IM_MAX(im)
+		    z2 = IM_MIN(im)
+		} else {
+		    len_stdline = SAMPLE_SIZE / NSAMPLE_LINES(cl)
+		    call zscale (im, z1, z2, CONTRAST(cl), SAMPLE_SIZE, 
+			len_stdline)
+		}
+	    }
+	}
+
+	if (strncmp (ZTRANS(cl), "auto", 1) == 0) {
+	    W_ZT(w1) = W_LINEAR
+	    # Calculate optimum z1 and z2 from image mode
+	    len_stdline = SAMPLE_SIZE / NSAMPLE_LINES(cl)
+	    call zscale (im, z1, z2, CONTRAST(cl), SAMPLE_SIZE, len_stdline)
+	    if (IM_PIXTYPE(im) == TY_SHORT) {
+	        if (short (z1) == short (z2))
+		    call error (0, 
+			"No range in data, ztrans=auto failed: z1 = z2")
+	    }
+	    else if (fp_equalr (z1, z2))
+		call error (0, "No range in data, ztrans=auto failed: z1=z2")
+	}
+
+	# Set the intensity extremes of the window descriptor
+	if (strncmp (ZTRANS(cl), "none", 1) == 0) {
+	    W_ZT(w1) = W_UNITARY
+	    W_ZS(w1) = IM_MIN(im)
+	    W_ZE(w1) = IM_MAX(im)
+	} else {
+	    W_ZS(w1) = z1
+	    W_ZE(w1) = z2
+	}
+end
diff --git a/pkg/plot/crtpict/sigl2.x b/pkg/plot/crtpict/sigl2.x
new file mode 100644
index 00000000..c1e1c9fc
--- /dev/null
+++ b/pkg/plot/crtpict/sigl2.x
@@ -0,0 +1,677 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<error.h>
+
+.help sigl2, sigl2_setup
+.nf ___________________________________________________________________________
+SIGL2 -- Get a line from a spatially scaled 2-dimensional image.  This procedure
+works like the regular IMIO get line procedure, but rescales the input
+2-dimensional image in either or both axes upon input.  If the magnification
+ratio required is greater than 0 and less than 2 then linear interpolation is
+used to resample the image.  If the magnification ratio is greater than or
+equal to 2 then the image is block averaged by the smallest factor which
+reduces the magnification to the range 0-2 and then interpolated back up to
+the desired size.  In some cases this will smooth the data slightly, but the
+operation is efficient and avoids aliasing effects.
+
+	si =	sigl2_setup (im, x1,x2,nx, y1,y2,ny)
+		sigl2_free (si)
+	ptr =	sigl2[sr] (si, linenumber)
+
+SIGL2_SETUP must be called to set up the transformations after mapping the
+image and before performing any scaled i/o to the image.  SIGL2_FREE must be
+called when finished to return buffer space.
+.endhelp ______________________________________________________________________
+
+# Scaled image descriptor for 2-dim images
+
+define	SI_LEN		15
+define	SI_MAXDIM	2		# images of 2 dimensions supported
+define	SI_NBUFS	3		# nbuffers used by SIGL2
+
+define	SI_IM		Memi[$1]	# pointer to input image header
+define	SI_GRID		Memi[$1+1+$2-1]	# pointer to array of X coords
+define	SI_NPIX		Memi[$1+3+$2-1]	# number of X coords
+define	SI_BAVG		Memi[$1+5+$2-1]	# X block averaging factor
+define	SI_INTERP	Memi[$1+7+$2-1]	# interpolate X axis
+define	SI_BUF		Memi[$1+9+$2-1]	# line buffers
+define	SI_TYBUF	Memi[$1+12]	# buffer type
+define	SI_XOFF		Memi[$1+13]	# offset in input image to first X
+define	SI_INIT		Memi[$1+14]	# YES until first i/o is done
+
+define	OUTBUF		SI_BUF($1,3)
+
+define	SI_TOL		(1E-5)		# close to a pixel
+define	INTVAL		(abs ($1 - nint($1)) < SI_TOL)
+define	SWAPI		{tempi=$2;$2=$1;$1=tempi}
+define	SWAPP		{tempp=$2;$2=$1;$1=tempp}
+define	NOTSET		(-9999)
+
+# SIGL2_SETUP -- Set up the spatial transformation for SIGL2[SR].  Compute
+# the block averaging factors (1 if no block averaging is required) and
+# the sampling grid points, i.e., pixel coordinates of the output pixels in
+# the input image.
+#
+# Valdes - Jan 9, 1985:
+#    Nx or ny can be 1 and blocking factors can be specified.
+
+pointer procedure sigl2_setup (im, px1, px2, nx, xblk, py1, py2, ny, yblk)
+
+pointer	im			# the input image
+real	px1, px2		# range in X to be sampled on an even grid
+int	nx			# number of output pixels in X
+int	xblk			# blocking factor in x
+real	py1, py2		# range in Y to be sampled on an even grid
+int	ny			# number of output pixels in Y
+int	yblk			# blocking factor in y
+
+int	npix, noldpix, nbavpix, i, j
+int	npts[SI_MAXDIM]		# number of output points for axis
+int	blksize[SI_MAXDIM]	# block averaging factor (npix per block)
+real	tau[SI_MAXDIM]		# tau = p(i+1) - p(i) in fractional pixels
+real	p1[SI_MAXDIM]		# starting pixel coords in each axis
+real	p2[SI_MAXDIM]		# ending pixel coords in each axis
+real	scalar, start
+pointer	si, gp
+
+begin
+	iferr (call calloc (si, SI_LEN, TY_STRUCT))
+	    call erract (EA_FATAL)
+
+	SI_IM(si) = im
+	SI_NPIX(si,1) = nx
+	SI_NPIX(si,2) = ny
+	SI_INIT(si) = YES
+
+	p1[1] = px1			# X = index 1
+	p2[1] = px2
+	npts[1] = nx
+	blksize[1] = xblk
+
+	p1[2] = py1			# Y = index 2
+	p2[2] = py2
+	npts[2] = ny
+	blksize[2] = yblk
+
+	# Compute block averaging factors if not defined.
+	# If there is only one pixel then the block average is the average
+	# between the first and last point.
+
+	do i = 1, SI_MAXDIM {
+	    if ((blksize[i] >= 1) && !IS_INDEFI(blksize[i])) {
+	        if (npts[i] == 1)
+		    tau[i] = 0.
+	        else
+	            tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+	    } else {
+		if (npts[i] == 1) {
+		    tau[i] = 0.
+		    blksize[i] = int (p2[i] - p1[i] + 1)
+	        } else {
+	            tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+	    	    if (tau[i] >= 2.0) {
+
+			# If nx or ny is not an integral multiple of the block
+			# averaging factor, noldpix is the next larger number
+			# which is an integral multiple.  When the image is
+			# block averaged pixels will be replicated as necessary
+			# to fill the last block out to this size.  
+
+			blksize[i] = int (tau[i])
+			npix = p2[i] - p1[i] + 1
+			noldpix = (npix+blksize[i]-1) / blksize[i] * blksize[i]
+			nbavpix = noldpix / blksize[i]
+			scalar = real (nbavpix - 1) / real (noldpix - 1)
+			p1[i] = (p1[i] - 1.0) * scalar + 1.0
+			p2[i] = (p2[i] - 1.0) * scalar + 1.0
+			tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+		    } else
+			blksize[i] = 1
+		}
+	    }
+	}
+
+	SI_BAVG(si,1) = blksize[1]
+	SI_BAVG(si,2) = blksize[2]
+
+	if (IS_INDEFI (xblk))
+	    xblk = blksize[1]
+	if (IS_INDEFI (yblk))
+	    yblk = blksize[2]
+
+	# Allocate and initialize the grid arrays, specifying the X and Y
+	# coordinates of each pixel in the output image, in units of pixels
+	# in the input (possibly block averaged) image.
+
+	do i = 1, SI_MAXDIM {
+	    # The X coordinate is special.  We do not want to read entire
+	    # input image lines if only a range of input X values are needed.
+	    # Since the X grid vector passed to ALUI (the interpolator) must
+	    # contain explicit offsets into the vector being interpolated,
+	    # we must generate interpolator grid points starting near 1.0.
+	    # The X origin, used to read the block averaged input line, is
+	    # given by XOFF.
+
+	    if (i == 1) {
+		SI_XOFF(si) = int (p1[i])
+		start = p1[1] - int (p1[i]) + 1.0
+	    } else
+		start = p1[i]
+
+	    # Do the axes need to be interpolated?
+	    if (INTVAL(start) && INTVAL(tau[i]))
+		SI_INTERP(si,i) = NO
+	    else
+		SI_INTERP(si,i) = YES
+
+	    # Allocate grid buffer and set the grid points.
+	    iferr (call malloc (gp, npts[i], TY_REAL))
+		call erract (EA_FATAL)
+	    SI_GRID(si,i) = gp
+	    do j = 0, npts[i]-1
+		Memr[gp+j] = start + (j * tau[i])
+	}
+
+	return (si)
+end
+
+
+# SIGL2_FREE -- Free storage associated with an image opened for scaled
+# input.  This does not close and unmap the image.
+
+procedure sigl2_free (si)
+
+pointer	si
+int	i
+
+begin
+	# Free SIGL2 buffers.
+	do i = 1, SI_NBUFS
+	    if (SI_BUF(si,i) != NULL)
+		call mfree (SI_BUF(si,i), SI_TYBUF(si))
+
+	# Free GRID buffers.
+	do i = 1, SI_MAXDIM
+	    if (SI_GRID(si,i) != NULL)
+		call mfree (SI_GRID(si,i), TY_REAL)
+
+	call mfree (si, TY_STRUCT)
+end
+
+
+# SIGL2S -- Get a line of type short from a scaled image.  Block averaging is
+# done by a subprocedure; this procedure gets a line from a possibly block
+# averaged image and if necessary interpolates it to the grid points of the
+# output line.
+
+pointer procedure sigl2s (si, lineno)
+
+pointer	si		# pointer to SI descriptor
+int	lineno
+
+pointer	rawline, tempp, gp
+int	i, buf_y[2], new_y[2], tempi, curbuf, altbuf
+int	npix, nblks_y, ybavg, x1, x2
+real	x, y, weight_1, weight_2
+pointer	si_blkavgs()
+errchk	si_blkavgs
+
+begin
+	npix = SI_NPIX(si,1)
+
+	# Determine the range of X (in pixels on the block averaged input image)
+	# required for the interpolator.
+
+	gp = SI_GRID(si,1)
+	x1 = SI_XOFF(si)
+	x = Memr[gp+npix-1]
+	x2 = x1 + int(x)
+	if (INTVAL(x))
+	    x2 = x2 - 1
+	x2 = max (x1 + 1, x2)
+
+	gp = SI_GRID(si,2)
+	y = Memr[gp+lineno-1]
+
+	# The following is an optimization provided for the case when it is
+	# not necessary to interpolate in either X or Y.  Block averaging is
+	# permitted.
+
+	if (SI_INTERP(si,1) == NO && SI_INTERP(si,2) == NO)
+	    return (si_blkavgs (SI_IM(si), x1, x2, int(y),
+		SI_BAVG(si,1), SI_BAVG(si,2)))
+
+	# If we are interpolating in Y two buffers are required, one for each
+	# of the two input image lines required to interpolate in Y.  The lines
+	# stored in these buffers are interpolated in X to the output grid but
+	# not in Y.  Both buffers are not required if we are not interpolating
+	# in Y, but we use them anyhow to simplify the code.
+
+	if (SI_INIT(si) == YES) {
+	    do i = 1, 2 {
+		if (SI_BUF(si,i) != NULL)
+		    call mfree (SI_BUF(si,i), SI_TYBUF(si))
+		call malloc (SI_BUF(si,i), npix, TY_SHORT)
+		SI_TYBUF(si) = TY_SHORT
+		buf_y[i] = NOTSET
+	    }
+	    if (OUTBUF(si) != NULL)
+		call mfree (OUTBUF(si), SI_TYBUF(si))
+	    call malloc (OUTBUF(si), npix, TY_SHORT)
+	    SI_INIT(si) = NO
+	}
+
+	# If the Y value of the new line is not in range of the contents of the
+	# current line buffers, refill one or both buffers.  To refill we must
+	# read a (possibly block averaged) input line and interpolate it onto
+	# the X grid.  The X and Y values herein are in the coordinate system
+	# of the (possibly block averaged) input image.
+
+	new_y[1] = int(y)
+	new_y[2] = int(y) + 1
+
+	# Get the pair of lines whose integral Y values form an interval
+	# containing the fractional Y value of the output line.  Sometimes the
+	# desired line will happen to be in the other buffer already, in which
+	# case we just have to swap buffers.  Often the new line will be the
+	# current line, in which case nothing is done.  This latter case occurs
+	# frequently when the magnification ratio is large.
+
+	curbuf = 1
+	altbuf = 2
+
+	do i = 1, 2 {
+	    if (new_y[i] == buf_y[i]) {
+		;
+	    } else if (new_y[i] == buf_y[altbuf]) {
+		SWAPP (SI_BUF(si,1), SI_BUF(si,2))
+		SWAPI (buf_y[1], buf_y[2])
+
+	    } else {
+		# Get line and interpolate onto output grid.  If interpolation
+		# is not required merely copy data out.  This code is set up
+		# to always use two buffers; in effect, there is one buffer of
+		# look ahead, even when Y[i] is integral.  This means that we
+		# will go out of bounds by one line at the top of the image.
+		# This is handled by copying the last line.
+
+		ybavg = SI_BAVG(si,2)
+		nblks_y = (IM_LEN (SI_IM(si), 2) + ybavg-1) / ybavg
+		if (new_y[i] <= nblks_y)
+		    rawline = si_blkavgs (SI_IM(si), x1, x2, new_y[i],
+			SI_BAVG(si,1), SI_BAVG(si,2))
+
+		if (SI_INTERP(si,1) == NO)
+		    call amovs (Mems[rawline], Mems[SI_BUF(si,i)], npix)
+		else {
+		    call aluis (Mems[rawline], Mems[SI_BUF(si,i)],
+			Memr[SI_GRID(si,1)], npix)
+		}
+
+		buf_y[i] = new_y[i]
+	    }
+
+	    SWAPI (altbuf, curbuf)
+	}
+
+	# We now have two line buffers straddling the output Y value,
+	# interpolated to the X grid of the output line.  To complete the
+	# bilinear interpolation operation we take a weighted sum of the two
+	# lines.  If the range from buf_y[1] to buf_y[2] is repeatedly
+	# interpolated in Y no additional i/o occurs and the linear
+	# interpolation operation (ALUI) does not have to be repeated (only the
+	# weighted sum is required).  If the distance of Y from one of the
+	# buffers is zero then we do not even have to take a weighted sum.
+	# This is not unusual because we may be called with a magnification
+	# of 1.0 in Y.
+
+	weight_1 = 1.0 - (y - buf_y[1])
+	weight_2 = 1.0 - weight_1
+
+	if (weight_2 < SI_TOL) 
+	    return (SI_BUF(si,1))
+	else if (weight_1 < SI_TOL)
+	    return (SI_BUF(si,2))
+	else {
+	    call awsus (Mems[SI_BUF(si,1)], Mems[SI_BUF(si,2)],
+		Mems[OUTBUF(si)], npix, weight_1, weight_2)
+	    return (OUTBUF(si))
+	}
+end
+
+
+# SI_BLKAVGS -- Get a line from a block averaged image of type short.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on.  If the length of an axis
+# is not an integral multiple of the block size then the last pixel in the
+# last block will be replicated to fill out the block; the average is still
+# defined even if a block is not full.
+
+pointer procedure si_blkavgs (im, x1, x2, y, xbavg, ybavg)
+
+pointer	im			# input image
+int	x1, x2			# range of x blocks to be read
+int	y			# y block to be read
+int	xbavg, ybavg		# X and Y block averaging factors
+
+short	snlines
+int	nblks_x, nblks_y, ncols, nlines, xoff, i, j
+int	first_line, nlines_in_sum, npix, nfull_blks, count
+real	sum
+pointer	sp, a, b
+pointer	imgs2s()
+errchk	imgs2s
+
+begin
+	call smark (sp)
+
+	ncols  = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xoff   = (x1 - 1) * xbavg + 1
+	npix   = min (ncols, xoff + (x2 - x1 + 1) * xbavg - 1)
+
+	if ((xbavg < 1) || (ybavg < 1))
+	    call error (1, "si_blkavg: illegal block size")
+	else if (x1 < 1 || x2 > ncols)
+	    call error (2, "si_blkavg: column index out of bounds")
+	else if ((xbavg == 1) && (ybavg == 1))
+	    return (imgs2s (im, xoff, xoff + npix - 1, y, y))
+
+	nblks_x = (npix   + xbavg-1) / xbavg
+	nblks_y = (nlines + ybavg-1) / ybavg
+
+	if (y < 1 || y > nblks_y)
+	    call error (2, "si_blkavg: block number out of range")
+
+	call salloc (b, nblks_x, TY_SHORT)
+
+	if (ybavg > 1) {
+	    call aclrs (Mems[b], nblks_x)
+	    nlines_in_sum = 0
+	}
+
+	# Read and accumulate all input lines in the block.
+	first_line = (y - 1) * ybavg + 1
+
+	do i = first_line, min (nlines, first_line + ybavg - 1) {
+	    # Get line from input image.
+	    a = imgs2s (im, xoff, xoff + npix - 1, i, i)
+
+	    # Block average line in X.
+	    if (xbavg > 1) {
+		# First block average only the full blocks.
+		nfull_blks = npix / xbavg
+		call abavs (Mems[a], Mems[a], nfull_blks, xbavg)
+
+		# Now average the final partial block, if any.
+		if (nfull_blks < nblks_x) {
+		    sum = 0.0
+		    count = 0
+		    do j = nfull_blks * xbavg + 1, npix {
+			sum = sum + Mems[a+j-1]
+			count = count + 1
+		    }
+		    Mems[a+nblks_x-1] = sum / count
+		}
+	    }
+
+	    # Add line into block sum.  Keep track of number of lines in sum
+	    # so that we can compute block average later.
+	    if (ybavg > 1) {
+		call aadds (Mems[a], Mems[b], Mems[b], nblks_x)
+		nlines_in_sum = nlines_in_sum + 1
+	    }
+	}
+
+	# Compute the block average in Y from the sum of all lines block
+	# averaged in X.  Overwrite buffer A, the buffer returned by IMIO.
+	# This is kosher because the block averaged line is never longer
+	# than an input line.
+
+	if (ybavg > 1) {
+	    snlines = short (nlines_in_sum)
+	    call adivks (Mems[b], snlines, Mems[a], nblks_x)
+	}
+
+	call sfree (sp)
+	return (a)
+end
+
+
+# SIGL2R -- Get a line of type real from a scaled image.  Block averaging is
+# done by a subprocedure; this procedure gets a line from a possibly block
+# averaged image and if necessary interpolates it to the grid points of the
+# output line.
+
+pointer procedure sigl2r (si, lineno)
+
+pointer	si		# pointer to SI descriptor
+int	lineno
+
+pointer	rawline, tempp, gp
+int	i, buf_y[2], new_y[2], tempi, curbuf, altbuf
+int	npix, nblks_y, ybavg, x1, x2
+real	x, y, weight_1, weight_2
+pointer	si_blkavgr()
+errchk	si_blkavgr
+
+begin
+	npix = SI_NPIX(si,1)
+
+	# Deterine the range of X (in pixels on the block averaged input image)
+	# required for the interpolator.
+
+	gp = SI_GRID(si,1)
+	x1 = SI_XOFF(si)
+	x = Memr[gp+npix-1]
+	x2 = x1 + int(x)
+	if (INTVAL(x))
+	    x2 = x2 - 1
+	x2 = max (x1 + 1, x2)
+
+	gp = SI_GRID(si,2)
+	y = Memr[gp+lineno-1]
+
+	# The following is an optimization provided for the case when it is
+	# not necessary to interpolate in either X or Y.  Block averaging is
+	# permitted.
+
+	if (SI_INTERP(si,1) == NO && SI_INTERP(si,2) == NO)
+	    return (si_blkavgr (SI_IM(si), x1, x2, int(y),
+		SI_BAVG(si,1), SI_BAVG(si,2)))
+
+	# If we are interpolating in Y two buffers are required, one for each
+	# of the two input image lines required to interpolate in Y.  The lines
+	# stored in these buffers are interpolated in X to the output grid but
+	# not in Y.  Both buffers are not required if we are not interpolating
+	# in Y, but we use them anyhow to simplify the code.
+
+	if (SI_INIT(si) == YES) {
+	    do i = 1, 2 {
+		if (SI_BUF(si,i) != NULL)
+		    call mfree (SI_BUF(si,i), SI_TYBUF(si))
+		call malloc (SI_BUF(si,i), npix, TY_REAL)
+		SI_TYBUF(si) = TY_REAL
+		buf_y[i] = NOTSET
+	    }
+	    if (OUTBUF(si) != NULL)
+		call mfree (OUTBUF(si), SI_TYBUF(si))
+	    call malloc (OUTBUF(si), npix, TY_REAL)
+	    SI_INIT(si) = NO
+	}
+
+	# If the Y value of the new line is not in range of the contents of the
+	# current line buffers, refill one or both buffers.  To refill we must
+	# read a (possibly block averaged) input line and interpolate it onto
+	# the X grid.  The X and Y values herein are in the coordinate system
+	# of the (possibly block averaged) input image.
+
+	new_y[1] = int(y)
+	new_y[2] = int(y) + 1
+
+	# Get the pair of lines whose integral Y values form an interval
+	# containing the fractional Y value of the output line.  Sometimes the
+	# desired line will happen to be in the other buffer already, in which
+	# case we just have to swap buffers.  Often the new line will be the
+	# current line, in which case nothing is done.  This latter case occurs
+	# frequently when the magnification ratio is large.
+
+	curbuf = 1
+	altbuf = 2
+
+	do i = 1, 2 {
+	    if (new_y[i] == buf_y[i]) {
+		;
+	    } else if (new_y[i] == buf_y[altbuf]) {
+		SWAPP (SI_BUF(si,1), SI_BUF(si,2))
+		SWAPI (buf_y[1], buf_y[2])
+
+	    } else {
+		# Get line and interpolate onto output grid.  If interpolation
+		# is not required merely copy data out.  This code is set up
+		# to always use two buffers; in effect, there is one buffer of
+		# look ahead, even when Y[i] is integral.  This means that we
+		# will go out of bounds by one line at the top of the image.
+		# This is handled by copying the last line.
+
+		ybavg = SI_BAVG(si,2)
+		nblks_y = (IM_LEN (SI_IM(si), 2) + ybavg-1) / ybavg
+		if (new_y[i] <= nblks_y)
+		    rawline = si_blkavgr (SI_IM(si), x1, x2, new_y[i],
+			SI_BAVG(si,1), SI_BAVG(si,2))
+
+		if (SI_INTERP(si,1) == NO)
+		    call amovr (Memr[rawline], Memr[SI_BUF(si,i)], npix)
+		else {
+		    call aluir (Memr[rawline], Memr[SI_BUF(si,i)],
+			Memr[SI_GRID(si,1)], npix)
+		}
+
+		buf_y[i] = new_y[i]
+	    }
+
+	    SWAPI (altbuf, curbuf)
+	}
+
+	# We now have two line buffers straddling the output Y value,
+	# interpolated to the X grid of the output line.  To complete the
+	# bilinear interpolation operation we take a weighted sum of the two
+	# lines.  If the range from buf_y[1] to buf_y[2] is repeatedly
+	# interpolated in Y no additional i/o occurs and the linear
+	# interpolation operation (ALUI) does not have to be repeated (only the
+	# weighted sum is required).  If the distance of Y from one of the
+	# buffers is zero then we do not even have to take a weighted sum.
+	# This is not unusual because we may be called with a magnification
+	# of 1.0 in Y.
+
+	weight_1 = 1.0 - (y - buf_y[1])
+	weight_2 = 1.0 - weight_1
+
+	if (weight_2 < SI_TOL) 
+	    return (SI_BUF(si,1))
+	else if (weight_1 < SI_TOL)
+	    return (SI_BUF(si,2))
+	else {
+	    call awsur (Memr[SI_BUF(si,1)], Memr[SI_BUF(si,2)],
+		Memr[OUTBUF(si)], npix, weight_1, weight_2)
+	    return (OUTBUF(si))
+	}
+end
+
+
+# SI_BLKAVGR -- Get a line from a block averaged image of type short.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on.  If the length of an axis
+# is not an integral multiple of the block size then the last pixel in the
+# last block will be replicated to fill out the block; the average is still
+# defined even if a block is not full.
+
+pointer procedure si_blkavgr (im, x1, x2, y, xbavg, ybavg)
+
+pointer	im			# input image
+int	x1, x2			# range of x blocks to be read
+int	y			# y block to be read
+int	xbavg, ybavg		# X and Y block averaging factors
+
+int	nblks_x, nblks_y, ncols, nlines, xoff, i, j
+int	first_line, nlines_in_sum, npix, nfull_blks, count
+real	sum
+pointer	sp, a, b
+pointer	imgs2r()
+errchk	imgs2r
+
+begin
+	call smark (sp)
+
+	ncols  = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xoff   = (x1 - 1) * xbavg + 1
+	npix   = min (ncols, xoff + (x2 - x1 + 1) * xbavg - 1)
+
+	if ((xbavg < 1) || (ybavg < 1))
+	    call error (1, "si_blkavg: illegal block size")
+	else if (x1 < 1 || x2 > ncols)
+	    call error (2, "si_blkavg: column index out of bounds")
+	else if ((xbavg == 1) && (ybavg == 1))
+	    return (imgs2r (im, xoff, xoff + npix - 1, y, y))
+
+	nblks_x = (npix   + xbavg-1) / xbavg
+	nblks_y = (nlines + ybavg-1) / ybavg
+
+	if (y < 1 || y > nblks_y)
+	    call error (2, "si_blkavg: block number out of range")
+
+	call salloc (b, nblks_x, TY_REAL)
+
+	if (ybavg > 1) {
+	    call aclrr (Memr[b], nblks_x)
+	    nlines_in_sum = 0
+	}
+
+	# Read and accumulate all input lines in the block.
+	first_line = (y - 1) * ybavg + 1
+
+	do i = first_line, min (nlines, first_line + ybavg - 1) {
+	    # Get line from input image.
+	    a = imgs2r (im, xoff, xoff + npix - 1, i, i)
+
+	    # Block average line in X.
+	    if (xbavg > 1) {
+		# First block average only the full blocks.
+		nfull_blks = npix / xbavg
+		call abavr (Memr[a], Memr[a], nfull_blks, xbavg)
+
+		# Now average the final partial block, if any.
+		if (nfull_blks < nblks_x) {
+		    sum = 0.0
+		    count = 0
+		    do j = nfull_blks * xbavg + 1, npix {
+			sum = sum + Memr[a+j-1]
+			count = count + 1
+		    }
+		    Memr[a+nblks_x-1] = sum / count
+		}
+	    }
+
+	    # Add line into block sum.  Keep track of number of lines in sum
+	    # so that we can compute block average later.
+	    if (ybavg > 1) {
+		call aaddr (Memr[a], Memr[b], Memr[b], nblks_x)
+		nlines_in_sum = nlines_in_sum + 1
+	    }
+	}
+
+	# Compute the block average in Y from the sum of all lines block
+	# averaged in X.  Overwrite buffer A, the buffer returned by IMIO.
+	# This is kosher because the block averaged line is never longer
+	# than an input line.
+
+	if (ybavg > 1)
+	    call adivkr (Memr[b], real(nlines_in_sum), Memr[a], nblks_x)
+
+	call sfree (sp)
+	return (a)
+end
diff --git a/pkg/plot/crtpict/t_crtpict.x b/pkg/plot/crtpict/t_crtpict.x
new file mode 100644
index 00000000..bd1887f9
--- /dev/null
+++ b/pkg/plot/crtpict/t_crtpict.x
@@ -0,0 +1,162 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include <error.h>
+include	<gset.h>
+include	<fset.h>
+include	<imhdr.h>
+include	"crtpict.h"
+
+# T_CRTPICT -- Code for the CRTPICT replacement.  Input to the procedure
+# is an IRAF image; output is a file of GKI metacode instructions, 
+# essentially x,y,z for each pixel on the dicomed.  The image intensities
+# are scaled to the dynamic range of the output device.  The image
+# is also scaled spatially, either expanded or reduced.  
+
+procedure t_crtpict ()
+
+bool	redir
+pointer	sp, cl, gp, im, command, image, word, title, output, ofile, dev
+int	cmd, stat, fd
+
+pointer	immap(), gopen()
+bool	clgetb(), streq()
+int	strncmp(), clgeti(), btoi(), fstati(), open(), getline()
+int	imtopenp(), list, imtgetim()
+real	clgetr()
+
+begin
+	call smark (sp)
+	call salloc (cl, LEN_CLPAR, TY_STRUCT)
+	call salloc (command, SZ_LINE, TY_CHAR)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (word, SZ_LINE, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (ofile, SZ_FNAME, TY_CHAR)
+	call salloc (dev, SZ_FNAME, TY_CHAR)
+
+	# If the input has been redirected, input is read from the named
+	# command file.  If not, each image name in the input template is
+	# plotted.
+	
+	call fseti (STDOUT, F_FLUSHNL, YES)
+	if (fstati (STDIN, F_REDIR) == YES) {
+	    call printf ("Input has been redirected\n")
+	    redir = true
+	    cmd = open ("STDIN", READ_ONLY, TEXT_FILE)
+	} else 
+	    list = imtopenp ("input")
+
+	# The user can "trap" the output metacode and intercept the
+	# spooling process if an output file is specified.
+	call clgstr ("output", Memc[output], SZ_FNAME)
+	if (!streq (Memc[output], "")) {
+	    call strcpy (Memc[output], Memc[ofile], SZ_FNAME)
+	    fd = open (Memc[ofile], NEW_FILE, BINARY_FILE)
+	} else
+	    fd = STDGRAPH
+
+	# Now get the parameters necessary to calculate z1, z2.
+	call clgstr ("ztrans", ZTRANS(cl), SZ_LINE)
+	if (strncmp (ZTRANS(cl), "auto", 1) == 0) {
+	    CONTRAST(cl) = clgetr ("contrast")
+	    NSAMPLE_LINES(cl) = clgeti ("nsample_lines")
+	} else if (strncmp (ZTRANS(cl), "min_max", 1) == 0) {
+	    Z1(cl) = clgetr ("z1")
+	    Z2(cl) = clgetr ("z2")
+	    if (abs (Z1(cl) - Z2(cl)) < EPSILON) {
+		CONTRAST(cl) = clgetr ("contrast")
+		if (abs (CONTRAST(cl)) - 1.0 > EPSILON)
+		    NSAMPLE_LINES(cl) = clgeti ("nsample_lines")
+	    }
+	} else if (strncmp (ZTRANS(cl), "user", 1) == 0)
+	    call clgstr ("lutfile", UFILE(cl), SZ_FNAME)
+
+	# Get parameters necessary to compute the spatial transformation.
+	FILL(cl) = btoi (clgetb ("auto_fill"))
+	if (FILL(cl) == NO) {
+	    XMAG(cl) = clgetr ("xmag")
+	    YMAG(cl) = clgetr ("ymag")
+	    if (XMAG(cl) < 0)
+		XMAG(cl) = 1 / abs (XMAG(cl))
+	    if (YMAG(cl) < 0)
+		YMAG(cl) = 1 / abs (YMAG(cl))
+	}
+
+	if (FILL(cl) == YES || XMAG(cl) - 1.0 > EPSILON || 
+	    YMAG(cl) - 1.0 > EPSILON) {
+	    # Find out how spatial scaling is to be done
+	    REPLICATE(cl) = btoi (clgetb ("replicate"))
+	    if (REPLICATE(cl) == NO) {
+	        # Get block averaging factors
+	        X_BA(cl) = clgetr ("x_blk_average")
+	        Y_BA(cl) = clgetr ("y_blk_average")
+	    }
+	} else
+	    REPLICATE(cl) = YES
+
+	# And for the plotting fractions:
+	PERIM(cl) = btoi (clgetb ("perimeter"))
+	GRAPHICS_FRACTION(cl) = clgetr ("graphics_fraction")
+	IMAGE_FRACTION(cl) = clgetr ("image_fraction")
+	GREYSCALE_FRACTION(cl) = clgetr ("greyscale_fraction")
+
+	# Get the output device name and determine the graphics pointer.
+	call clgstr ("device", Memc[dev], SZ_FNAME)
+	gp = gopen (Memc[dev], NEW_FILE, fd)
+
+	# Loop over commands until EOF
+	repeat {
+	    if (redir) {
+		if (getline (STDIN, Memc[command]) == EOF)
+		    break
+		call sscan (Memc[command])
+		    call gargwrd (Memc[word], SZ_LINE)
+		if (!streq (Memc[word], "plot")) {
+		    # Pixel window has been stored as WCS 2
+		    call gseti (gp, G_WCS, 2)
+		    call gscan (gp, Memc[command])
+		    next
+		} else 
+		    call gargwrd (Memc[image], SZ_FNAME)
+	    } else {
+		stat = imtgetim (list, Memc[image], SZ_FNAME)
+		if (stat == EOF)
+		    break
+	    }
+
+	    # Open the input image; if an error occurs, go to next image in list
+	    iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    iferr (call crt_plot_image (gp, im, Memc[image], cl))
+		call erract (EA_WARN)
+
+	    # Add title to metacode file
+	    call sprintf (Memc[title], SZ_LINE, "CRTPICT: %s")
+		call pargstr (IM_TITLE(im))
+	    call gmftitle (gp, Memc[title])
+	    call imunmap (im)
+
+	    # Clear frame for next picture (unless plotting to terminal)
+	    if (strncmp (Memc[dev], "stdgraph", 4) == 0)
+	        call gflush (gp)
+	    else
+	        call gclear (gp)
+	}
+
+	# Clean up and close files
+	call gclose (gp)
+	call close (fd)
+
+	if (redir)
+	    call close (cmd)
+	else
+	    call imtclose (list)
+
+	call sfree (sp)
+	call flush (STDOUT)
+end
diff --git a/pkg/plot/crtpict/tweakndc.x b/pkg/plot/crtpict/tweakndc.x
new file mode 100644
index 00000000..1a8f674b
--- /dev/null
+++ b/pkg/plot/crtpict/tweakndc.x
@@ -0,0 +1,66 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# CRT_TWEAK_NDC -- Alter the input ndc endpoints so that the ratio
+# of device_pixels to image_pixels is an integer.  This procedure insures
+# an integer replication (or decimation) factor.  For replication, 
+# ndevice_pixels_spanned / nimage_values_output = an_integer.  For
+# decimation, ndevice_pixels_spanned / nimage_values_output = 1 / an_integer.
+# The NDC coordinates returned also represent integer device pixels.
+
+procedure crt_tweak_ndc (nvalues, ndc_start, ndc_end, device_res)
+
+int	nvalues			# Number of image values to output
+real	ndc_start, ndc_end	# NDC endpoints - changed on output
+int	device_res		# Full device resolution
+
+int	ndevice_elements, first_device_element, last_device_element, int_extra
+int	dev_pixel_1, dev_pixel_2, desired_dev_elements, desired_inverse
+real	ndc_extra, extra, real_extra
+double	gki_tweak
+
+begin
+	gki_tweak = double (32767) / double (32768)
+	if (int (ndc_end) == 1)
+	    ndc_end = gki_tweak
+	first_device_element = ndc_start * device_res 
+	last_device_element  = ndc_end * device_res
+	ndevice_elements = (last_device_element - first_device_element) + 1
+
+	if (mod (ndevice_elements, nvalues) != 0) {
+	    # Calculate amount to be altered
+	    real_extra = real (ndevice_elements) / real (nvalues)
+	    if (real_extra > 1.0) {
+	        # Tweak to get an integer replication factor
+	        int_extra  = ndevice_elements / nvalues
+	        extra = real ((real_extra - int_extra) * nvalues)
+	        ndc_extra = extra / device_res
+	        ndc_start = ndc_start + (ndc_extra / 2.0)
+	        ndc_end   = ndc_end   - (ndc_extra / 2.0)
+	    } else {
+		# Tweak to get an integer decimation factor
+		real_extra = real (nvalues) / real (ndevice_elements)
+		desired_inverse = int (real_extra) + 1
+		desired_dev_elements = nvalues / desired_inverse
+		extra = desired_dev_elements - ndevice_elements
+		ndc_extra = real (extra) / real (device_res)
+		ndc_start = ndc_start - (ndc_extra / 2.0)
+		ndc_end   = ndc_end   + (ndc_extra / 2.0)
+	    }
+	}
+
+	# Now have ndc coordinates of starting and ending pixel such
+	# that the replication or decimation factor is
+	# an integer.  Now insure that the ndc coordinates refer to 
+	# integer device pixels  so that truncation later in the
+	# processing doesn't alter this replication factor.  In what
+	# follows, note that dev_pixel_1 and dev_pixel_2
+	# are 0-based; dev_pixel_1 is the first pixel to be filled, and 
+	# dev_pixel_2 is the first pixel NOT to be filled, in accordance
+	# with Richard's notes.
+
+	dev_pixel_1 = ndc_start * device_res
+	dev_pixel_2 = (ndc_end  * device_res) + 1
+
+	ndc_start = real (dev_pixel_1) / real (device_res) / gki_tweak
+	ndc_end = real (dev_pixel_2) / real (device_res) / gki_tweak
+end
diff --git a/pkg/plot/crtpict/wdes.h b/pkg/plot/crtpict/wdes.h
new file mode 100644
index 00000000..60d89f21
--- /dev/null
+++ b/pkg/plot/crtpict/wdes.h
@@ -0,0 +1,33 @@
+# Window descriptor structure.
+
+define	LEN_WDES	(5+(W_MAXWC+1)*LEN_WC+40)
+define	LEN_WC		9			# 4=[XbXeYbYe]+2=tr_type[xy]
+define	W_MAXWC		5			# max world coord systems
+define	W_SZIMSECT	79			# image section string
+
+define	W_DEVICE	Memi[$1]
+define	W_FRAME		Memi[$1+1]		# device frame number
+define	W_XRES		Memi[$1+2]		# device resolution, x
+define	W_YRES		Memi[$1+3]		# device resolution, y
+define	W_WC		($1+$2*LEN_WC+5)	# ptr to coord descriptor
+define	W_IMSECT	Memc[P2C($1+50)]
+
+# Fields of the WC coordinate descriptor, a substructure of the window
+# descriptor.  "W_XB(W_WC(w,0))" is the XB field of wc 0 of window W.
+
+define	W_XS		Memr[P2R($1)]		# starting X value
+define	W_XE		Memr[P2R($1+1)]		# ending X value
+define	W_XT		Memi[$1+2]		# X transformation type
+define	W_YS		Memr[P2R($1+3)]		# starting Y value
+define	W_YE		Memr[P2R($1+4)]		# ending Y value
+define	W_YT		Memi[$1+5]		# Y transformation type
+define	W_ZS		Memr[P2R($1+6)]		# starting Z value (greyscale)
+define	W_ZE		Memr[P2R($1+7)]		# ending Z value
+define	W_ZT		Memi[$1+8]		# Z transformation type
+
+# Types of coordinate and greyscale transformations.
+
+define	W_UNITARY	0			# values map without change
+define	W_LINEAR	1			# linear mapping
+define	W_LOG		2			# logarithmic mapping
+define	W_USER		3			# user supplied lut values
diff --git a/pkg/plot/crtpict/xformimage.x b/pkg/plot/crtpict/xformimage.x
new file mode 100644
index 00000000..0400cc91
--- /dev/null
+++ b/pkg/plot/crtpict/xformimage.x
@@ -0,0 +1,117 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<mach.h>
+include	<gset.h>
+include	"wdes.h"
+include	"crtpict.h"
+
+# CRT_TRANSFORM_IMAGE -- Map an image into the output device.  In general this
+# involves independent linear transformations in the X, Y, and Z (greyscale)
+# dimensions.  If a spatial dimension is larger than the display window then
+# the image is block averaged.  If a spatial dimension or a block averaged
+# dimension is smaller than the display window then linear interpolation is
+# used to expand the image.  The input image is accessed via IMIO; the 
+# transformed image is output to a greyscale device via GIO.    
+
+procedure crt_transform_image (gp, im, wdes, cl)
+
+pointer	gp			# graphics descriptor for output
+pointer	im			# input image
+pointer	wdes			# graphics window descriptor
+pointer	cl
+
+real	ndc_xs,ndc_xe,ndc_ys,ndc_ye	# NDC of output device window
+real	px1, px2, py1, py2	# image coords in fractional image pixels
+real	pxsize, pysize		# size of image section in fractional pixels
+real	wxcenter, wycenter	# center of device window in frac device pixels
+real	xmag, ymag		# x,y magnification ratios
+pointer	w0, w1			# world coord systems 0 (NDC) and 1 (pixel)
+int	ndev_cols, ndev_rows, nx_output, ny_output
+int	ggeti()
+errchk	ggeti, gseti, gsview, gswind, draw_perimeter
+errchk	crt_map_image, crt_replicate_image
+
+begin
+	# Compute pointers to WCS 0 and 1.
+	w0 = W_WC(wdes,0)
+	w1 = W_WC(wdes,1)
+
+	call printf ("%s: %s\n")
+	    call pargstr (W_IMSECT(wdes))
+	    call pargstr (IM_TITLE(im))
+
+	ndev_cols = ggeti (gp, "xr")
+	ndev_rows = ggeti (gp, "yr")
+
+	# Compute X and Y magnification ratios required to map image into
+	# the device window.
+
+	xmag = ((W_XE(w0) - W_XS(w0)) * ndev_cols) / (W_XE(w1) - W_XS(w1))
+	ymag = ((W_YE(w0) - W_YS(w0)) * ndev_rows) / (W_YE(w1) - W_YS(w1))
+
+	# Compute the coordinates of the image section to be displayed.
+	# This is not necessarily the same as WCS 1 since the WCS coords
+	# need not be inbounds, but they are integral pixels.
+
+	px1 = max (1.0,		        real (int (W_XS(w1) + 0.5)))
+	px2 = min (real (IM_LEN(im,1)), real (int (W_XE(w1) + 0.5)))
+	py1 = max (1.0,		        real (int (W_YS(w1) + 0.5)))
+	py2 = min (real (IM_LEN(im,2)), real (int (W_YE(w1) + 0.5)))
+
+	# The extent of the output image in NDC coords
+	pxsize = ((px2 - px1) + 1.0) / ndev_cols
+	pysize = ((py2 - py1) + 1.0) / ndev_rows
+
+	# The NDC coordinates that map to the central image pixel output
+	wxcenter = (W_XE(w0) + W_XS(w0)) / 2.0
+	wycenter = (W_YE(w0) + W_YS(w0)) / 2.0
+
+	# Now compute the NDC coordinates of the output device that the
+	# image will occupy.  
+	ndc_xs = max (W_XS(w0), wxcenter - (pxsize / 2.0 * xmag))
+	ndc_xe = max (ndc_xs, min (W_XE(w0), ndc_xs + (pxsize * xmag)))
+	ndc_ys = max (W_YS(w0), wycenter - (pysize / 2.0 * ymag))
+	ndc_ye = max (ndc_ys, min (W_YE(w0), ndc_ys + (pysize * ymag)))
+
+	# To avoid possible truncation errors down the line, make sure
+	# the ndc coordinates passed to the output procedures represent
+	# integer pixels.
+	ndc_xs = real (int (ndc_xs * ndev_cols)) / ndev_cols
+	ndc_xe = real (int (ndc_xe * ndev_cols)) / ndev_cols
+	ndc_ys = real (int (ndc_ys * ndev_rows)) / ndev_rows
+	ndc_ye = real (int (ndc_ye * ndev_rows)) / ndev_rows
+
+	# Output the image data in WCS 0.  The number of image pixels that
+	# will be put out across the device is calculated first.  
+
+	call gseti (gp, G_WCS, 0)
+	if (REPLICATE(cl) == YES) {
+	    nx_output = (int(px2) - int(px1)) + 1
+	    ny_output = (int(py2) - int(py1)) + 1
+	} else {
+	    # Image pixels will be scaled to number of device pixels 
+	    nx_output = ((ndc_xe * ndev_cols) - (ndc_xs * ndev_cols)) + 1
+	    ny_output = ((ndc_ye * ndev_rows) - (ndc_ys * ndev_rows)) + 1
+	    nx_output = nx_output / X_BA(cl)
+	    ny_output = ny_output / Y_BA(cl)
+	}
+
+	# Tweak the ndc coordinates to insure integer replication factors.
+	# This may change the ndc coordinates.
+	call crt_tweak_ndc (nx_output, ndc_xs, ndc_xe, ndev_cols)
+	call crt_tweak_ndc (ny_output, ndc_ys, ndc_ye, ndev_rows)
+
+	# Call routine that actually puts out the pixels row by row
+	call crt_map_image (im, gp, px1,px2,py1,py2, ndc_xs,ndc_xe,
+	       ndc_ys,ndc_ye, nx_output, ny_output, W_ZS(w1),W_ZE(w1),W_ZT(w1),
+	           cl)
+
+	# Change to pixel coordinates and draw perimeter axes if requested
+	call gseti  (gp, G_WCS, 2)
+	call gsview (gp, ndc_xs, ndc_xe, ndc_ys, ndc_ye)
+	call gswind (gp, px1, px2, py1, py2)
+
+	if (PERIM(cl) == YES)
+	    call draw_perimeter (gp)
+end
diff --git a/pkg/plot/crtpict/xyscale.x b/pkg/plot/crtpict/xyscale.x
new file mode 100644
index 00000000..094b78ee
--- /dev/null
+++ b/pkg/plot/crtpict/xyscale.x
@@ -0,0 +1,90 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include <imhdr.h>
+include	<error.h>
+include	"wdes.h"
+include	"crtpict.h"
+
+# CRT_XY_SCALE -- Calculate the spatial transformation parameters and store
+# them in the "wdes" window descriptor.  These are used to map the image
+# to the graphics device.
+
+procedure crt_xy_scale (gp, cl, im, wdes)
+
+pointer	gp
+pointer	cl
+pointer	im
+pointer	wdes
+
+int	w1, w0, ndev_rows, ndev_cols
+real	xcenter, ycenter, xsize, ysize, pxsize, pysize, xscale, yscale
+real	ystart, yend
+int	ggeti()
+errchk	ggeti
+
+begin
+	ndev_rows = ggeti (gp, "yr")
+	ndev_cols = ggeti (gp, "xr")
+
+	# Determine the maximum display window available for mapping the
+	# image.  This is a function of the user specified "fractions";
+	# the outer 6% of the user specified image area is used for annotation.
+
+	xcenter = (CRT_XE + CRT_XS) / 2.0
+	xsize   = CRT_XE - CRT_XS
+	ystart  = CRT_YS + (CRT_YE - CRT_YS) * GRAPHICS_FRACTION(cl)
+	yend    = ystart + (CRT_YE - CRT_YS) * IMAGE_FRACTION(cl)
+	ycenter = (yend + ystart) / 2.0
+	ysize   = yend - ystart
+
+	# Set device window limits in normalized device coordinates.
+	# World coord system 0 is used for the device window.
+
+	w0 = W_WC(wdes,0)
+	W_XS(w0) = xcenter - (xsize / 2.0)
+	W_XE(w0) = xcenter + (xsize / 2.0)
+	W_YS(w0) = ycenter - (ysize / 2.0)
+	W_YE(w0) = ycenter + (ysize / 2.0)
+
+	# Determine X and Y scaling ratios required to map the image into the
+	# normalized display window.  
+
+	if (FILL(cl) == YES) {
+	    # Compute scale in units of NDC viewport coords per image pixel.
+	    xscale = xsize / max (1, (IM_LEN(im,1) - 1))
+	    yscale = ysize / max (1, (IM_LEN(im,2) - 1))
+
+	    # Image scaled by same factor in x and y with FILL = YES.
+	    if (xscale < yscale)
+		yscale = xscale
+	    else
+		xscale = yscale
+
+	} else {
+	    # From the user specified magnification factors, compute the scale
+	    # in units of NDC coords per "effective" display pixels.
+
+	    xscale = 1.0 / ((ndev_cols - 1) / XMAG(cl))
+	    yscale = 1.0 / ((ndev_rows - 1) / YMAG(cl))
+	}
+
+	# Determine the image pixels that map to the available device
+	# viewport.  The starting and ending pixels calculated from
+	# xscale and yscale may be in the image interior or reference 
+	# beyond the bounds of the image.  These endpoints are tuned
+	# further in procedure transform_image.
+
+	w1 = W_WC(wdes,1)
+	pxsize = xsize / xscale
+	pysize = ysize / yscale
+
+	W_XS(w1) = (IM_LEN(im,1) - 1) / 2.0 + 1 - (pxsize / 2.0)
+	W_XE(w1) = W_XS(w1) + pxsize
+	W_YS(w1) = (IM_LEN(im,2) - 1) / 2.0 + 1 - (pysize / 2.0)
+	W_YE(w1) = W_YS(w1) + pysize
+
+	# All spatial transformations are linear.
+	W_XT(w1) = W_LINEAR
+	W_YT(w1) = W_LINEAR
+end
diff --git a/pkg/plot/crtpict/zscale.x b/pkg/plot/crtpict/zscale.x
new file mode 100644
index 00000000..aa0b925d
--- /dev/null
+++ b/pkg/plot/crtpict/zscale.x
@@ -0,0 +1,441 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+.help zscale
+.nf ___________________________________________________________________________
+ZSCALE -- Compute the optimal Z1, Z2 (range of greyscale values to be
+displayed) of an image.  For efficiency a statistical subsample of an image
+is used.  The pixel sample evenly subsamples the image in x and y.  The entire
+image is used if the number of pixels in the image is smaller than the desired
+sample.
+
+The sample is accumulated in a buffer and sorted by greyscale value.
+The median value is the central value of the sorted array.  The slope of a
+straight line fitted to the sorted sample is a measure of the standard
+deviation of the sample about the median value.  Our algorithm is to sort
+the sample and perform an iterative fit of a straight line to the sample,
+using pixel rejection to omit gross deviants near the endpoints.  The fitted
+straight line is the transfer function used to map image Z into display Z.
+If more than half the pixels are rejected the full range is used.  The slope
+of the fitted line is divided by the user-supplied contrast factor and the
+final Z1 and Z2 are computed, taking the origin of the fitted line at the
+median value.
+.endhelp ______________________________________________________________________
+
+define	MIN_NPIXELS	5		# smallest permissible sample
+define	MAX_REJECT	0.5		# max frac. of pixels to be rejected
+define	GOOD_PIXEL	0		# use pixel in fit
+define	BAD_PIXEL	1		# ignore pixel in all computations
+define	REJECT_PIXEL	2		# reject pixel after a bit
+define	KREJ		2.5		# k-sigma pixel rejection factor
+define	MAX_ITERATIONS	5		# maximum number of fitline iterations
+
+
+# ZSCALE -- Sample the image and compute Z1 and Z2.
+
+procedure zscale (im, z1, z2, contrast, optimal_sample_size, len_stdline)
+
+pointer	im			# image to be sampled
+real	z1, z2			# output min and max greyscale values
+real	contrast		# adj. to slope of transfer function
+int	optimal_sample_size	# desired number of pixels in sample
+int	len_stdline		# optimal number of pixels per line
+
+int	npix, minpix, ngoodpix, center_pixel, ngrow
+real	zmin, zmax, median, temp
+real	zstart, zslope
+pointer	sample, left
+int	zsc_sample_image(), zsc_fit_line()
+
+begin
+	# Subsample the image.
+	npix = zsc_sample_image (im, sample, optimal_sample_size, len_stdline)
+	center_pixel = max (1, (npix + 1) / 2)
+
+	# Sort the sample, compute the minimum, maximum, and median pixel
+	# values.
+
+	call asrtr (Memr[sample], Memr[sample], npix)
+	zmin = Memr[sample]
+	zmax = Memr[sample+npix-1]
+
+	# The median value is the average of the two central values if there 
+	# are an even number of pixels in the sample.
+
+	left = sample + center_pixel - 1
+	if (mod (npix, 2) == 1 || center_pixel >= npix)
+	    median = Memr[left]
+	else
+	    median = (Memr[left] + Memr[left+1]) / 2
+
+	# Fit a line to the sorted sample vector.  If more than half of the
+	# pixels in the sample are rejected give up and return the full range.
+	# If the user-supplied contrast factor is not 1.0 adjust the scale
+	# accordingly and compute Z1 and Z2, the y intercepts at indices 1 and
+	# npix.
+
+	minpix = max (MIN_NPIXELS, int (npix * MAX_REJECT))
+	ngrow = max (1, nint (npix * .01))
+	ngoodpix = zsc_fit_line (Memr[sample], npix, zstart, zslope,
+	    KREJ, ngrow, MAX_ITERATIONS)
+
+	if (ngoodpix < minpix) {
+	    z1 = zmin
+	    z2 = zmax
+	} else {
+	    zslope = zslope / abs (contrast)
+	    z1 = max (zmin, median - (center_pixel - 1) * zslope)
+	    z2 = min (zmax, median + (npix - center_pixel) * zslope)
+	    if (contrast < 0) {
+		# Negative contrast desired, flip z1 and z2
+		temp = z1
+		z1 = z2
+		z2 = temp
+	    }
+	}
+
+	call mfree (sample, TY_REAL)
+end
+
+
+# ZSC_SAMPLE_IMAGE -- Extract an evenly gridded subsample of the pixels from
+# a two-dimensional image into a one-dimensional vector.
+
+int procedure zsc_sample_image (im, sample, optimal_sample_size, len_stdline)
+
+pointer	im			# image to be sampled
+pointer	sample			# output vector containing the sample
+int	optimal_sample_size	# desired number of pixels in sample
+int	len_stdline		# optimal number of pixels per line
+
+int	ncols, nlines, col_step, line_step, maxpix, line
+int	opt_npix_per_line, npix_per_line
+int	opt_nlines_in_sample, min_nlines_in_sample, max_nlines_in_sample
+pointer	op
+pointer	imgl2r()
+
+begin
+	ncols  = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+
+	# Compute the number of pixels each line will contribute to the sample,
+	# and the subsampling step size for a line.  The sampling grid must
+	# span the whole line on a uniform grid.
+
+	opt_npix_per_line = min (ncols, len_stdline)
+	col_step = (ncols + opt_npix_per_line-1) / opt_npix_per_line
+	npix_per_line = (ncols + col_step-1) / col_step
+
+	# Compute the number of lines to sample and the spacing between lines.
+	# We must ensure that the image is adequately sampled despite its
+	# size, hence there is a lower limit on the number of lines in the
+	# sample.  We also want to minimize the number of lines accessed when
+	# accessing a large image, because each disk seek and read is expensive.
+	# The number of lines extracted will be roughly the sample size divided
+	# by len_stdline, possibly more if the lines are very short.
+
+	min_nlines_in_sample = max (1, optimal_sample_size / len_stdline)
+	opt_nlines_in_sample = max(min_nlines_in_sample, min(nlines,
+	    (optimal_sample_size + npix_per_line-1) / npix_per_line))
+	line_step = max (1, nlines / (opt_nlines_in_sample))
+	max_nlines_in_sample = (nlines + line_step-1) / line_step
+
+	# Allocate space for the output vector.  Buffer must be freed by our
+	# caller.
+
+	maxpix = npix_per_line * max_nlines_in_sample
+	call malloc (sample, maxpix, TY_REAL)
+
+	# Extract the vector.
+	op = sample
+	do line = (line_step + 1) / 2, nlines, line_step {
+	    call zsc_subsample (Memr[imgl2r(im,line)], Memr[op],
+		npix_per_line, col_step)
+	    op = op + npix_per_line
+	    if (op - sample + npix_per_line > maxpix)
+		break
+	}
+
+	return (op - sample)
+end
+
+
+# ZSC_SUBSAMPLE -- Subsample an image line.  Extract the first pixel and
+# every "step"th pixel thereafter for a total of npix pixels.
+
+procedure zsc_subsample (a, b, npix, step)
+
+real	a[ARB]
+real	b[npix]
+int	npix, step
+int	ip, i
+
+begin
+	if (step <= 1)
+	    call amovr (a, b, npix)
+	else {
+	    ip = 1
+	    do i = 1, npix {
+		b[i] = a[ip]
+		ip = ip + step
+	    }
+	}
+end
+
+
+# ZSC_FIT_LINE -- Fit a straight line to a data array of type real.  This is
+# an iterative fitting algorithm, wherein points further than ksigma from the
+# current fit are excluded from the next fit.  Convergence occurs when the
+# next iteration does not decrease the number of pixels in the fit, or when
+# there are no pixels left.  The number of pixels left after pixel rejection
+# is returned as the function value.
+
+int procedure zsc_fit_line (data, npix, zstart, zslope, krej, ngrow, maxiter)
+
+real	data[npix]		# data to be fitted
+int	npix			# number of pixels before rejection
+real	zstart			# Z-value of pixel data[1]	(output)
+real	zslope			# dz/pixel			(output)
+real	krej			# k-sigma pixel rejection factor
+int	ngrow			# number of pixels of growing
+int	maxiter			# max iterations
+
+int	i, ngoodpix, last_ngoodpix, minpix, niter
+real	xscale, z0, dz, x, z, mean, sigma, threshold
+double	sumxsqr, sumxz, sumz, sumx, rowrat
+pointer	sp, flat, badpix, normx
+int	zsc_reject_pixels(), zsc_compute_sigma()
+
+begin
+	call smark (sp)
+
+	if (npix <= 0)
+	    return (0)
+	else if (npix == 1) {
+	    zstart = data[1]
+	    zslope = 0.0
+	    return (1)
+	} else
+	    xscale = 2.0 / (npix - 1)
+
+	# Allocate a buffer for data minus fitted curve, another for the
+	# normalized X values, and another to flag rejected pixels.
+
+	call salloc (flat, npix, TY_REAL)
+	call salloc (normx, npix, TY_REAL)
+	call salloc (badpix, npix, TY_SHORT)
+	call aclrs (Mems[badpix], npix)
+
+	# Compute normalized X vector.  The data X values [1:npix] are
+	# normalized to the range [-1:1].  This diagonalizes the lsq matrix
+	# and reduces its condition number.
+
+	do i = 0, npix - 1
+	    Memr[normx+i] = i * xscale - 1.0
+
+	# Fit a line with no pixel rejection.  Accumulate the elements of the
+	# matrix and data vector.  The matrix M is diagonal with
+	# M[1,1] = sum x**2 and M[2,2] = ngoodpix.  The data vector is
+	# DV[1] = sum (data[i] * x[i]) and DV[2] = sum (data[i]).
+
+	sumxsqr = 0
+	sumxz = 0
+	sumx = 0
+	sumz = 0
+
+	do i = 1, npix {
+	    x = Memr[normx+i-1]
+	    z = data[i]
+	    sumxsqr = sumxsqr + (x ** 2)
+	    sumxz   = sumxz + z * x
+	    sumz    = sumz + z
+	}
+	# Solve for the coefficients of the fitted line.
+	z0 = sumz / npix
+	dz = sumxz / sumxsqr
+
+	# Iterate, fitting a new line in each iteration.  Compute the flattened
+	# data vector and the sigma of the flat vector.  Compute the lower and
+	# upper k-sigma pixel rejection thresholds.  Run down the flat array
+	# and detect pixels to be rejected from the fit.  Reject pixels from
+	# the fit by subtracting their contributions from the matrix sums and
+	# marking the pixel as rejected.
+
+	ngoodpix = npix
+	minpix = max (MIN_NPIXELS, int (npix * MAX_REJECT))
+
+	for (niter=1;  niter <= maxiter;  niter=niter+1) {
+	    last_ngoodpix = ngoodpix
+
+	    # Subtract the fitted line from the data array.
+	    call zsc_flatten_data (data, Memr[flat], Memr[normx], npix, z0, dz)
+
+	    # Compute the k-sigma rejection threshold.  In principle this
+	    # could be more efficiently computed using the matrix sums
+	    # accumulated when the line was fitted, but there are problems with
+	    # numerical stability with that approach.
+
+	    ngoodpix = zsc_compute_sigma (Memr[flat], Mems[badpix], npix,
+		mean, sigma)
+	    threshold = sigma * krej
+
+	    # Detect and reject pixels further than ksigma from the fitted
+	    # line.
+	    ngoodpix = zsc_reject_pixels (data, Memr[flat], Memr[normx],
+		Mems[badpix], npix, sumxsqr, sumxz, sumx, sumz, threshold,
+		ngrow)
+
+	    # Solve for the coefficients of the fitted line.  Note that after
+	    # pixel rejection the sum of the X values need no longer be zero.
+
+	    if (ngoodpix > 0) {
+		rowrat = sumx / sumxsqr
+		z0 = (sumz - rowrat * sumxz) / (ngoodpix - rowrat * sumx)
+		dz = (sumxz - z0 * sumx) / sumxsqr
+	    }
+
+	    if (ngoodpix >= last_ngoodpix || ngoodpix < minpix)
+		break
+	}
+
+	# Transform the line coefficients back to the X range [1:npix].
+	zstart = z0 - dz
+	zslope = dz * xscale
+
+	call sfree (sp)
+	return (ngoodpix)
+end
+
+
+# ZSC_FLATTEN_DATA -- Compute and subtract the fitted line from the data array,
+# returned the flattened data in FLAT.
+
+procedure zsc_flatten_data (data, flat, x, npix, z0, dz)
+
+real	data[npix]		# raw data array
+real	flat[npix]		# flattened data  (output)
+real	x[npix]			# x value of each pixel
+int	npix			# number of pixels
+real	z0, dz			# z-intercept, dz/dx of fitted line
+int	i
+
+begin
+	do i = 1, npix
+	    flat[i] = data[i] - (x[i] * dz + z0)
+end
+
+
+# ZSC_COMPUTE_SIGMA -- Compute the root mean square deviation from the
+# mean of a flattened array.  Ignore rejected pixels.
+
+int procedure zsc_compute_sigma (a, badpix, npix, mean, sigma)
+
+real	a[npix]			# flattened data array
+short	badpix[npix]		# bad pixel flags (!= 0 if bad pixel)
+int	npix
+real	mean, sigma		# (output)
+
+real	pixval
+int	i, ngoodpix
+double	sum, sumsq, temp
+
+begin
+	sum = 0
+	sumsq = 0
+	ngoodpix = 0
+
+	# Accumulate sum and sum of squares.
+	do i = 1, npix
+	    if (badpix[i] == GOOD_PIXEL) {
+		pixval = a[i]
+		ngoodpix = ngoodpix + 1
+		sum = sum + pixval
+		sumsq = sumsq + pixval ** 2
+	    }
+
+	# Compute mean and sigma.
+	switch (ngoodpix) {
+	case 0:
+	    mean = INDEF
+	    sigma = INDEF
+	case 1:
+	    mean = sum
+	    sigma = INDEF
+	default:
+	    mean = sum / ngoodpix
+	    temp = sumsq / (ngoodpix - 1) - sum**2 / (ngoodpix * (ngoodpix - 1))
+	    if (temp < 0)		# possible with roundoff error
+		sigma = 0.0
+	    else
+		sigma = sqrt (temp)
+	}
+
+	return (ngoodpix)
+end
+
+
+# ZSC_REJECT_PIXELS -- Detect and reject pixels more than "threshold" greyscale
+# units from the fitted line.  The residuals about the fitted line are given
+# by the "flat" array, while the raw data is in "data".  Each time a pixel
+# is rejected subtract its contributions from the matrix sums and flag the
+# pixel as rejected.  When a pixel is rejected reject its neighbors out to
+# a specified radius as well.  This speeds up convergence considerably and
+# produces a more stringent rejection criteria which takes advantage of the
+# fact that bad pixels tend to be clumped.  The number of pixels left in the
+# fit is returned as the function value.
+
+int procedure zsc_reject_pixels (data, flat, normx, badpix, npix,
+				 sumxsqr, sumxz, sumx, sumz, threshold, ngrow)
+
+real	data[npix]		# raw data array
+real	flat[npix]		# flattened data array
+real	normx[npix]		# normalized x values of pixels
+short	badpix[npix]		# bad pixel flags (!= 0 if bad pixel)
+int	npix
+double	sumxsqr,sumxz,sumx,sumz	# matrix sums
+real	threshold		# threshold for pixel rejection
+int	ngrow			# number of pixels of growing
+
+int	ngoodpix, i, j
+real	residual, lcut, hcut
+double	x, z
+
+begin
+	ngoodpix = npix
+	lcut = -threshold
+	hcut = threshold
+
+	do i = 1, npix
+	    if (badpix[i] == BAD_PIXEL)
+		ngoodpix = ngoodpix - 1
+	    else {
+		residual = flat[i]
+		if (residual < lcut || residual > hcut) {
+		    # Reject the pixel and its neighbors out to the growing
+		    # radius.  We must be careful how we do this to avoid
+		    # directional effects.  Do not turn off thresholding on
+		    # pixels in the forward direction; mark them for rejection
+		    # but do not reject until they have been thresholded.
+		    # If this is not done growing will not be symmetric.
+
+		    do j = max(1,i-ngrow), min(npix,i+ngrow) {
+			if (badpix[j] != BAD_PIXEL) {
+			    if (j <= i) {
+				x = normx[j]
+				z = data[j]
+				sumxsqr = sumxsqr - (x ** 2)
+				sumxz = sumxz - z * x
+				sumx = sumx - x
+				sumz = sumz - z
+				badpix[j] = BAD_PIXEL
+				ngoodpix = ngoodpix - 1
+			    } else
+				badpix[j] = REJECT_PIXEL
+			}
+		    }
+		}
+	    }
+
+	return (ngoodpix)
+end