Repatch (from linux) of OSX IRAF

16 files changed, 1719 insertions, 0 deletions

diff --git a/sys/gio/glabax/README b/sys/gio/glabax/README
new file mode 100644
index 00000000..4c9f9ad5
--- /dev/null
+++ b/sys/gio/glabax/README
@@ -0,0 +1 @@
+GLABAX -- GIO axis drawing and labelling package.
diff --git a/sys/gio/glabax/glabax.h b/sys/gio/glabax/glabax.h
new file mode 100644
index 00000000..070918ec
--- /dev/null
+++ b/sys/gio/glabax/glabax.h
@@ -0,0 +1,46 @@
+# GLABAX.H -- Axis drawing and labelling.
+
+define	SZ_FORMAT		19
+define	SZ_LABEL		19
+define	MAX_LINEARITY		1.0			# no log scaling if gt
+define	LEFT_BORDER		9			# nchars at l|r edge
+define	BOTTOM_BORDER		5			# nlines at bottom edge
+define	Y_LABELOFFSET		7			# Y label dist from axis
+define	MAX_SZTITLEBLOCK	0.5			# max sztitleblock, NDC
+define	MIN_NTITLELINES		2			# min lines in titleblk
+define	TOL			(EPSILONR*10.0)
+
+define	LEN_AX			85
+define	AX_POS			Memd[P2D($1)+$2-1]	# tick coords
+define	AX_DRAWME		Memi[$1+4]		# draw this axis
+define	AX_HORIZONTAL		Memi[$1+5]		# axis is horizontal
+define	AX_SCALING		Memi[$1+6]		# type of scaling
+define	AX_DRAWTICKS		Memi[$1+7]		# draw the ticks
+define	AX_START		Memr[P2R($1+8+$2-1)]	# axis starts here
+define	AX_END			Memr[P2R($1+10+$2-1)]	# axis ends here
+define	AX_TICK1		Memr[P2R($1+12+$2-1)]	# first tick is here
+define	AX_STEP			Memr[P2R($1+14+$2-1)]	# offset between ticks
+define	AX_ISTEP		Memr[P2R($1+16+$2-1)]	# intial offset
+define	AX_KSTEP		Memr[P2R($1+18)]	# step scalar at majors
+define	AX_IKSTEP		Memr[P2R($1+19)]	# initial kstep
+define	AX_NMINOR		Memi[$1+20]		# nminor ticks
+define	AX_NLEFT		Memi[$1+21]		# nminor to next major
+define	AX_INLEFT		Memi[$1+22]		# initial nleft
+define	AX_NDIGITS		Memi[$1+23]		# ndigits of precision
+define	AX_MINORTICK		Memr[P2R($1+24+$2-1)]	# offset to draw minor
+define	AX_MAJORTICK		Memr[P2R($1+26+$2-1)]	# offset to draw major
+define	AX_MINORWIDTH		Memr[P2R($1+28)]	# minor tick linewidth
+define	AX_MAJORWIDTH		Memr[P2R($1+29)]	# major tick linewidth
+define	AX_LABELTICKS		Memi[$1+30]		# draw tick labels
+define	AX_TICKLABELOFFSET	Memr[P2R($1+31+$2-1)]	# offset to ticklabel
+define	AX_TICKLABELSIZE	Memr[P2R($1+33)]	# char size of ticklabel
+define	AX_TICKLABELCOLOR	Memi[$1+34]		# char size of ticklabel
+define	AX_TICKCOLOR		Memi[$1+35]		# grid between ticks
+define	AX_AXISLABELSIZE	Memr[P2R($1+36)]	# char size axislabel
+define	AX_AXISLABELCOLOR	Memi[$1+37]		# char size axislabel
+define	AX_AXISWIDTH		Memr[P2R($1+38)]	# axis linewidth
+define	AX_AXISCOLOR		Memi[$1+39]		# axis linewidth
+define	AX_GRIDCOLOR		Memi[$1+40]		# grid between ticks
+
+define	AX_TICKLABELPOS		Memc[P2C($1+45)]	# gtext format
+define	AX_TICKFORMAT		Memc[P2C($1+65)]	# numeric format
diff --git a/sys/gio/glabax/glabax.x b/sys/gio/glabax/glabax.x
new file mode 100644
index 00000000..0c30021b
--- /dev/null
+++ b/sys/gio/glabax/glabax.x
@@ -0,0 +1,264 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gset.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLABAX -- Draw and label the axes of the plot (normally the viewport
+# boundary).  This is done in two steps.  First we compute all the required
+# parameters, and then we draw and label the axes.  Up to four axes can be
+# drawn.  To simplify matters, all four axes are treated equally and
+# independently.  The axes are drawn a tick at a time in world coordinates.
+
+procedure glabax (gp, title, xlabel, ylabel)
+
+pointer	gp			# graphics descriptor
+char	title[ARB]		# plot title (may be more than one line)
+char	xlabel[ARB]		# X axis label
+char	ylabel[ARB]		# Y axis label
+
+char	label[SZ_LABEL]
+int	axis, wcs, ntitlelines, ip, major_tick
+int	save_plcolor, save_txcolor, save_facolor
+int	save_pltype, save_clip, save_txfont
+real	xv[4], yv[4], x1, x2, y1, y2
+real	save_plwidth, save_txsize
+real	dx, dy, x, y, sx, sy, scalar, wc, wstep
+pointer	sp, axes[4], ax, w
+
+real	gstatr()
+bool	ttygetb()
+int	gstati(), glb_gettick()
+errchk	glb_setup, gadraw, grdraw, gamove, gtext
+errchk	glb_label_axis, glb_plot_title, glb_gettick
+
+begin
+	call smark (sp)
+	call salloc (axes[1], LEN_AX, TY_STRUCT)
+	call salloc (axes[2], LEN_AX, TY_STRUCT)
+	call salloc (axes[3], LEN_AX, TY_STRUCT)
+	call salloc (axes[4], LEN_AX, TY_STRUCT)
+
+	wcs = GP_WCS(gp)
+	w   = GP_WCSPTR(gp,wcs)
+
+	# Count the number of lines in the title block.
+	ntitlelines = 0
+	if (title[1] != EOS) {
+	    for (ip=1;  title[ip] != EOS;  ip=ip+1)
+		if (title[ip] == '\n' && title[ip+1] != EOS)
+		    ntitlelines = ntitlelines + 1
+	    ntitlelines = ntitlelines + 1
+	}
+	ntitlelines = max (ntitlelines, GP_NTITLELINES(gp))
+
+	# Fix the coordinates systems and set the axis drawing parameters.
+	# The number of lines in the title block is needed to determine how
+	# much space to allow at the top of the screen.
+
+	call glb_setup (gp, axes, ntitlelines, xlabel, ylabel)
+
+	# Save the values of any user parameters we must change while drawing
+	# the axes.
+
+	save_pltype  = gstati (gp, G_PLTYPE)
+	save_plwidth = gstatr (gp, G_PLWIDTH)
+	save_plcolor = gstati (gp, G_PLCOLOR)
+	save_txfont  = gstati (gp, G_TXFONT)
+	save_txsize  = gstatr (gp, G_TXSIZE)
+	save_txcolor = gstati (gp, G_TXCOLOR)
+	save_facolor = gstati (gp, G_FACOLOR)
+	save_clip    = WCS_CLIP(w)
+
+	# Prepare the background.
+	if (ttygetb (GP_TTY(gp), "fa") &&
+	    GP_FRAMECOLOR(gp) != 0 && GP_FRAMEDRAWN(gp) == NO) {
+
+	    call ggview (gp, x1, x2, y1, y2)
+	    call gseti (gp, G_WCS, 0)
+	    call gseti (gp, G_CLIP, NO)
+
+	    xv[1] = 0.0;  yv[1] = 0.0
+	    xv[2] = 1.0;  yv[2] = 0.0
+	    xv[3] = 1.0;  yv[3] = 1.0
+	    xv[4] = 0.0;  yv[4] = 1.0
+	    call gseti (gp, G_FACOLOR, GP_FRAMECOLOR(gp))
+	    call gfill (gp, xv, yv, 4, GF_SOLID)
+
+	    xv[1] = x1;  yv[1] = y1
+	    xv[2] = x2;  yv[2] = y1
+	    xv[3] = x2;  yv[3] = y2
+	    xv[4] = x1;  yv[4] = y2
+	    call gseti (gp, G_FACOLOR, 0)
+	    call gfill (gp, xv, yv, 4, GF_SOLID)
+
+	    call gseti (gp, G_CLIP, save_clip)
+	    call gseti (gp, G_WCS, wcs)
+	    GP_FRAMEDRAWN(gp) = YES
+	}
+
+	# Draw and label the four axes.  First set the linetype and linewidth
+	# to be used to draw the axes and ticks; these may be different than
+	# that used to plot the data.  Draws are preferred to moves to minimize
+	# the number of polylines needed to draw the axis.  An axis is drawn
+	# by moving to the start of the axis, drawing each tick in sequence,
+	# and then moving to the end of the axis.  Tick labels are drawn at
+	# the major ticks if required.  The axes and ticks must be drawn in
+	# world coords to get the proper scaling.  Clipping is turned off while
+	# drawing the axes to avoid clipping portions of the axes due to small
+	# floating point errors.
+
+	call gseti (gp, G_PLTYPE, 1)
+	call gseti (gp, G_CLIP, NO)
+	call gseti (gp, G_TXFONT, GT_BOLD)
+
+	do axis = 1, 4 {
+	    ax = axes[axis]
+	    if (AX_DRAWME(ax) == NO)
+		next
+
+#  call eprintf ("axis %d: tick1=(%g,%g) istep=(%g,%g) kstep=%g\n")
+#  call pargi (axis)
+#  call pargr (AX_TICK1(ax,1)); call pargr (AX_TICK1(ax,2))
+#  call pargr (AX_ISTEP(ax,1)); call pargr (AX_ISTEP(ax,2))
+#  call pargr (AX_IKSTEP(ax))
+#  call eprintf ("\tstart=(%g,%g) end=(%g,%g)\n")
+#  call pargr (AX_START(ax,1)); call pargr (AX_START(ax,2))
+#  call pargr (AX_END(ax,1)); call pargr (AX_END(ax,2))
+#  call eprintf ("nminor=%d, inleft=%d, minortick=(%g,%g), majortick=(%g,%g)\n")
+#  call pargi (AX_NMINOR(ax)); call pargi (AX_INLEFT(ax))
+#  call pargr (AX_MINORTICK(ax,1)); call pargr (AX_MINORTICK(ax,2))
+#  call pargr (AX_MAJORTICK(ax,1)); call pargr (AX_MAJORTICK(ax,2))
+
+	    # Set the axis linewidth and move to the start of the axis.
+	    call gsetr (gp, G_PLWIDTH, AX_AXISWIDTH(ax))
+	    call gseti (gp, G_PLCOLOR, AX_AXISCOLOR(ax))
+	    call gamove (gp, AX_START(ax,1), AX_START(ax,2))
+
+	    # Draw the axis and label the major ticks if so indicated.
+	    # First set flag to initialize glb_gettick.
+
+	    AX_NLEFT(ax) = -1
+	    while (glb_gettick (gp, ax, x, y, major_tick) != EOF) {
+
+		# Advance to the next tick.
+		call gsetr (gp, G_PLWIDTH, AX_AXISWIDTH(ax))
+		call gseti (gp, G_PLCOLOR, AX_AXISCOLOR(ax))
+		call gadraw (gp, x, y)
+
+		if (major_tick == YES) {
+		    # Draw a major tick.
+
+		    call gsetr (gp, G_PLWIDTH, AX_MAJORWIDTH(ax))
+		    call gseti (gp, G_PLCOLOR, AX_TICKCOLOR(ax))
+		    dx = AX_MAJORTICK(ax,1)
+		    dy = AX_MAJORTICK(ax,2)
+		    call grdraw (gp, dx, dy)
+		    call grdraw (gp, -dx, -dy)
+
+		    if (AX_LABELTICKS(ax) == YES) {
+			# Get the tick label position in NDC coords.  World
+			# coords cannot be used for an offset outside the
+			# viewport as the coords might be indefinite if log
+			# scaling.
+
+			call gseti (gp, G_WCS, 0)
+			call gcurpos (gp, sx, sy)
+			dx = AX_TICKLABELOFFSET(ax,1)
+			dy = AX_TICKLABELOFFSET(ax,2)
+
+			# Format the numeric tick label string.  The scalar
+			# multiplier is used to compute the step size between
+			# major ticks.
+
+			scalar = AX_NMINOR(ax) + 1.0
+			if (AX_HORIZONTAL(ax) == YES) {
+			    wc = x
+			    wstep = AX_STEP(ax,1) * scalar
+			} else {
+			    wc = y
+			    wstep = AX_STEP(ax,2) * scalar
+			}
+
+			# Draw the label string.
+			call gsetr (gp, G_TXSIZE, AX_TICKLABELSIZE(ax))
+			call gseti (gp, G_TXCOLOR, AX_TICKLABELCOLOR(ax))
+
+			# If log scaling, label the ticks in log units.
+			if (AX_SCALING(ax) == LINEAR) {
+			    call glb_encode (wc, label, SZ_LABEL,
+				AX_TICKFORMAT(ax), wstep)
+			    call gtext (gp, sx + dx, sy + dy, label,
+				AX_TICKLABELPOS(ax))
+			} else {
+			    call glb_loglab (gp, sx+dx, sy+dy, wc,
+				AX_TICKLABELPOS(ax), AX_SCALING(ax))
+			}
+
+			# Leave the pen back at the base of the tick.
+			call gamove (gp, sx, sy)
+			call gseti (gp, G_WCS, wcs)
+		    } 
+
+		} else {
+		    # Draw a minor tick.
+
+		    dx = AX_MINORTICK(ax,1)
+		    dy = AX_MINORTICK(ax,2)
+
+		    call gsetr (gp, G_PLWIDTH, AX_MINORWIDTH(ax))
+		    call gseti (gp, G_PLCOLOR, AX_TICKCOLOR(ax))
+		    call grdraw (gp, dx, dy)
+		    call grdraw (gp, -dx, -dy)
+		}
+	    }
+
+	    # Draw line segment from last tick to the end of the axis.
+	    call gadraw (gp, AX_END(ax,1), AX_END(ax,2))
+
+	    # Flush the graphics output.  When working interactively, this
+	    # gives the user something to watch while we generate the rest
+	    # of the plot.
+
+	    if (AX_NMINOR(ax) > 0)
+		call gflush (gp)
+	}
+
+	# Draw grid between major ticks.
+	if (GL_DRAWGRID (GP_XAP(gp)) == YES) {
+	    call gseti (gp, G_PLCOLOR, AX_GRIDCOLOR(axes[3]))
+	    call glb_drawgrid (gp, axes[3], axes[2])
+	}
+	if (GL_DRAWGRID (GP_YAP(gp)) == YES) {
+	    call gseti (gp, G_PLCOLOR, AX_GRIDCOLOR(axes[1]))
+	    call glb_drawgrid (gp, axes[1], axes[4])
+	}
+
+	# Label the X and Y axes.
+	do axis = 1, 4 {
+	    ax = axes[axis]
+	    if (AX_DRAWME(ax) == YES && AX_LABELTICKS(ax) == YES) {
+		call gseti (gp, G_TXCOLOR, AX_AXISLABELCOLOR(ax))
+		call glb_label_axis (gp, ax, xlabel, ylabel)
+	    }
+	}
+
+	# Draw plot title block.
+	call gseti (gp, G_TXCOLOR, GP_TITLECOLOR(gp))
+	call glb_plot_title (gp, title, ntitlelines)
+
+	# Restore the parameters we were originally called with.
+	call gseti (gp, G_WCS, wcs)
+	call gseti (gp, G_CLIP,    save_clip)
+	call gseti (gp, G_PLTYPE,  save_pltype)
+	call gsetr (gp, G_PLWIDTH, save_plwidth)
+	call gseti (gp, G_PLCOLOR, save_plcolor)
+	call gsetr (gp, G_TXSIZE,  save_txsize)
+	call gseti (gp, G_TXFONT,  save_txfont)
+	call gseti (gp, G_TXCOLOR, save_txcolor)
+	call gseti (gp, G_FACOLOR, save_facolor)
+
+	call gflush (gp)
+	call sfree (sp)
+end
diff --git a/sys/gio/glabax/glbencode.x b/sys/gio/glabax/glbencode.x
new file mode 100644
index 00000000..cbed6875
--- /dev/null
+++ b/sys/gio/glabax/glbencode.x
@@ -0,0 +1,66 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	"glabax.h"
+
+# GLB_ENCODE -- Encode a floating point number as a character string for a
+# tick label.  We have to be careful how we do this, since on the one hand
+# we want the most concise label possible (e.g., 500 not 500.00) but on the
+# other we must provide enough precision to discriminate between ticks that
+# are close together (e.g., 500.02 and 500.04).  The extra information is
+# given by the "ndigits" argument, which was calculated knowing the range
+# and step at setup time.
+
+procedure glb_encode (x, out, maxch, format, step)
+
+real	x			# number to be encoded
+char	out[ARB]		# output string
+int	maxch			# max chars out
+char	format[ARB]		# sprintf format
+real	step			# tick spacing
+
+int	ip, op
+real	nicex
+define	trim_ 91
+
+begin
+	# Test for the zero tick, to avoid tick labels that look like the
+	# machine epsilon.
+
+	if (abs (x / step) < TOL)
+	    nicex = 0
+	else
+	    nicex = x
+
+	# Encode number.
+	call sprintf (out, maxch, format)
+	    call pargr (nicex)
+
+	# Lop off any insignificant trailing zeros or periods.  Watch out for
+	# trailing zeros in exponential format, e.g., "1.0E10".
+
+	for (ip=1;  out[ip] != EOS;  ip=ip+1)
+	    if (out[ip] == 'E' || out[ip] == 'D')
+		goto trim_
+
+	for (ip=ip-1;  ip > 1 && out[ip] == '0';  ip=ip-1)
+	    ;
+	if (ip > 1 && out[ip] == '.')
+	    ip = ip - 1
+	if (ip >= 1)
+	    out[ip+1] = EOS
+
+	# Lop off any insignificant leading zeros, but be sure to leave at
+	# least one digit.
+trim_
+	for (op=1;  out[op] == '-' || out[op] == '+';  op=op+1)
+	    ;
+	for (ip=op;  out[ip] == '0' && out[ip+1] != EOS;  ip=ip+1)
+	    ;
+	while (out[ip] != EOS) {
+	    out[op] = out[ip]
+	    op = op + 1
+	    ip = ip + 1
+	}
+	out[op] = EOS
+end
diff --git a/sys/gio/glabax/glbfind.x b/sys/gio/glabax/glbfind.x
new file mode 100644
index 00000000..b9ff3975
--- /dev/null
+++ b/sys/gio/glabax/glbfind.x
@@ -0,0 +1,339 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_FIND_TICKS -- Find the optimal positions for the tick marks on an axis.
+# If rounding is enabled, extend the WCS out to the next tick outside the
+# boundary on either end of the axis.  Since this routine may modify the WCS
+# it must be called before any other routines (e.g., glb_setaxes).  Our task
+# is to position the major ticks in world coordinates at round numbers, e.g.,
+# 10, 20, 30 for a linear scale or 10, 100, 1000 for a log scale.  The minor
+# ticks are evenly distributed over the range between the major ticks.  If
+# log scaling is in use the step size between ticks will change by an order
+# of magnitude (by the factor KSTEP) in each decade of the log scale, i.e.,
+# at each major tick.  All tick positions and offsets are output in world
+# coordinates.
+
+procedure glb_find_ticks (gp, ap, ax1, ax2, angle)
+
+pointer	gp			# graphics descriptor
+pointer	ap			# axis parameters (from graphics descriptor)
+pointer	ax1, ax2		# axis descriptors (output)
+int	angle			# axis orientation, 0 or 90 degrees
+
+pointer	w
+int	logflag, nminor, scaling, t1, t2
+real	char_height, char_width, wctick, tval
+real	p1, p2, tp1, tp2, wcp1, tick1, step, minor_step, length
+
+bool	fp_equalr()
+int	gt_ndigits()
+real	ggetr(), elogr(), aelogr(), glb_minorstep()
+
+begin
+	w = GP_WCSPTR (gp, GP_WCS(gp))
+
+	# Start by zeroing the AX structures so that we do not have to zero
+	# fields explicitly.  This is a bit tricky because we are implicitly
+	# setting fields that are not named, but it saves time and space.
+
+	call aclri (Memi[ax1], LEN_AX)
+	call aclri (Memi[ax2], LEN_AX)
+
+	# If ticks are not to be drawn or if there are fewer than 2 ticks then
+	# we are done.
+
+	if (GL_NMAJOR(ap) <= 2)
+	    return
+
+	# Call the tick placement algorithm to determine where to put the major
+	# tick marks.  The output of this block are the variables P1 and P1,
+	# the world coords of the ends of the axis, and TICK1 and STEP, the
+	# world coords of the first tick and separation in world coords between
+	# major ticks.
+
+	if (angle == 0) {
+	    p1 = WCS_WX1(w)
+	    p2 = WCS_WX2(w)
+	    scaling = WCS_XTRAN(w)
+	} else {
+	    p1 = WCS_WY1(w)
+	    p2 = WCS_WY2(w)
+	    scaling = WCS_YTRAN(w)
+	}
+
+	AX_SCALING(ax1) = scaling
+	AX_SCALING(ax2) = scaling
+
+	if (scaling == LOG) {
+	    p1 = log10 (p1)
+	    p2 = log10 (p2)
+	    logflag = YES
+	} else if (scaling == ELOG) {
+	    p1 = elogr (p1)
+	    p2 = elogr (p2)
+	    logflag = YES
+	} else
+	    logflag = NO
+
+	# Call the tick placement algorithm.
+	call gtickr (p1, p2, GL_NMAJOR(ap), logflag, tick1, step)
+
+	# If rounding is enabled, extend the WCS out to the next major tick
+	# position outward on either end.  Always round if log scaling.
+
+	if (GL_ROUND(ap) == YES || scaling != LINEAR) {
+	    if (!fp_equalr (p1, tick1)) {
+		tick1 = tick1 - step
+		p1 = tick1
+	    }
+
+	    length = (p2 - p1) / step
+	    if (!fp_equalr (p1 + int(length) * step, p2))
+		p2 = p1 + (int(length) + 1) * step
+
+	    if (scaling == ELOG) {
+		tp1 = aelogr (p1)
+		tp2 = aelogr (p2)
+	    } else if (scaling == LOG) {
+		tp1 = 10.0 ** p1
+		tp2 = 10.0 ** p2
+	    } else {
+		tp1 = p1
+		tp2 = p2
+	    }
+
+	    if (angle == 0) {
+		WCS_WX1(w) = tp1
+		WCS_WX2(w) = tp2
+	    } else {
+		WCS_WY1(w) = tp1
+		WCS_WY2(w) = tp2
+	    }
+
+	    GP_WCSSTATE(gp) = MODIFIED
+	}
+
+	# Compute the coords of the axis endpoint and of the first tick in world
+	# coords.
+
+	if (scaling == LINEAR) {
+	    wctick = tick1
+	    wcp1 = p1
+	} else if (scaling == LOG) {
+	    wctick = 10.0 ** tick1
+	    wcp1 = 10.0 ** p1
+	} else {
+	    wctick = aelogr (tick1)
+	    wcp1 = aelogr (p1)
+	}
+
+	# Compute the number of minor ticks.  If we are log scaling there
+	# are either no minor ticks or 8 minor ticks.  If the scaling is
+	# linear the tick placement algorithm is used to compute the best
+	# number of minor ticks, using GL_NMINOR as a close estimate.  If
+	# NMINOR is negative automatic tick selection is disabled and exactly
+	# abs(NMINOR) ticks will be drawn.  If NMINOR is zero no minor ticks
+	# are drawn.
+
+	if (GL_NMINOR(ap) == 0)			# no minor ticks
+	    nminor = 0
+	else if (logflag == YES)		# log scaling
+	    nminor = 8
+	else {
+	    minor_step = glb_minorstep (tick1, tick1+step, GL_NMINOR(ap))
+	    nminor = nint (abs (step / minor_step)) - 1
+	}
+
+	AX_NMINOR(ax1) = nminor
+	AX_NMINOR(ax2) = nminor
+
+	# Compute the step size in world coords between minor ticks and the
+	# number of minor ticks to be drawn initially until the first major
+	# tick (tick1) is reached.  Note that for ELOG scaling the minor
+	# step size and number of minor ticks are different in the range
+	# +-10 (which is linear) than elsewhere, but we ignore that here.
+
+	if (scaling == LINEAR) {
+	    minor_step = step / (nminor + 1)
+	    AX_INLEFT(ax1) = abs (int ((wctick - wcp1) / minor_step))
+	} else {
+	    t1 = nint (tick1)
+	    t2 = nint (tick1 + step)
+	    if (scaling == LOG)
+		minor_step = (10.0 ** t2 - 10.0 ** t1) / 9.
+	    else
+		minor_step = (aelogr(real(t2)) - aelogr(real(t1))) / 9.
+	    if (nminor == 0)
+		minor_step = minor_step * 9.
+	    AX_INLEFT(ax1) = 0
+	}
+
+	AX_INLEFT(ax2) = AX_INLEFT(ax1)
+
+	# Set KSTEP, the adjustment to the step size at each major tick.  This
+	# is always 1.0 if the scale is linear.  Set KSTEP to negative if ELOG
+	# scaling, to tell the drawing code to invert kstep (.1->10 or 10->.1)
+	# when passing through the origin (necessary for ELOG scaling).  The
+	# sign is not otherwise significant.  If heading toward the origin
+	# initially then KSTEP is inverted for ELOG scaling vs LOG scaling.
+
+	if (scaling == LINEAR) {
+	    AX_IKSTEP(ax1) = 1.0
+	} else if (scaling == ELOG) {
+	    tval = p1
+	    if (abs (tval + step) > abs(t1))
+		AX_IKSTEP(ax1) = -10.0
+	    else
+		AX_IKSTEP(ax1) = -0.1
+	} else
+	    AX_IKSTEP(ax1) = 10.0 ** step
+	AX_IKSTEP(ax2) = AX_IKSTEP(ax1)
+
+	# Set those parameters which differ depending on whether the axis is
+	# horizontal or vertical.
+
+	if (angle == 0) {
+	    AX_TICK1(ax1,1) = wctick - (AX_INLEFT(ax1) * minor_step)
+	    AX_TICK1(ax2,1) = wctick - (AX_INLEFT(ax2) * minor_step)
+
+	    if (GL_SETAXISPOS(ap) == YES) {
+		AX_TICK1(ax1,2) = GL_AXISPOS1(ap)
+		AX_TICK1(ax2,2) = GL_AXISPOS2(ap)
+	    } else {
+		AX_TICK1(ax1,2) = WCS_WY1(w)
+		AX_TICK1(ax2,2) = WCS_WY2(w)
+	    }
+
+	    AX_ISTEP(ax1,1) = minor_step
+	    AX_ISTEP(ax2,1) = minor_step
+
+	    char_height = ggetr (gp, "ch")
+	    if (char_height < EPSILON)
+		char_height = DEF_CHARHEIGHT
+	    char_height = char_height * GL_TICKLABELSIZE(ap)
+
+	    AX_TICKLABELOFFSET(ax2,2) = 0.5 * char_height
+	    AX_TICKLABELOFFSET(ax1,2) = -AX_TICKLABELOFFSET(ax2,2)
+
+	    # Set gtext format for tick labels.
+	    call strcpy ("hj=c,vj=t", AX_TICKLABELPOS(ax1), SZ_FORMAT)
+	    call strcpy ("hj=c,vj=b", AX_TICKLABELPOS(ax2), SZ_FORMAT)
+
+	} else {
+	    if (GL_SETAXISPOS(ap) == YES) {
+		AX_TICK1(ax1,1) = GL_AXISPOS1(ap)
+		AX_TICK1(ax2,1) = GL_AXISPOS2(ap)
+	    } else {
+		AX_TICK1(ax1,1) = WCS_WX1(w)
+		AX_TICK1(ax2,1) = WCS_WX2(w)
+	    }
+
+	    AX_TICK1(ax1,2) = wctick - (AX_INLEFT(ax1) * minor_step)
+	    AX_TICK1(ax2,2) = wctick - (AX_INLEFT(ax2) * minor_step)
+
+	    AX_ISTEP(ax1,2) = minor_step
+	    AX_ISTEP(ax2,2) = minor_step
+
+	    char_width = ggetr (gp, "cw")
+	    if (char_width < EPSILON)
+		char_width = DEF_CHARWIDTH
+	    char_width = char_width * GL_TICKLABELSIZE(ap)
+
+	    AX_TICKLABELOFFSET(ax2,1) = 0.5 * char_width
+	    AX_TICKLABELOFFSET(ax1,1) = -AX_TICKLABELOFFSET(ax2,1)
+
+	    call strcpy ("hj=r,vj=c", AX_TICKLABELPOS(ax1), SZ_FORMAT)
+	    call strcpy ("hj=l,vj=c", AX_TICKLABELPOS(ax2), SZ_FORMAT)
+	}
+
+	# Set the tick parameters that are identical for the two axes and
+	# which do not depend on whether the axis is horizontal or vertical.
+
+	AX_DRAWTICKS(ax1)	= GL_DRAWTICKS(ap)
+	AX_DRAWTICKS(ax2)	= GL_DRAWTICKS(ap)
+	AX_TICKLABELSIZE(ax1)	= GL_TICKLABELSIZE(ap)
+	AX_TICKLABELSIZE(ax2)	= GL_TICKLABELSIZE(ap)
+	AX_TICKLABELCOLOR(ax1)	= GL_TICKLABELCOLOR(ap)
+	AX_TICKLABELCOLOR(ax2)	= GL_TICKLABELCOLOR(ap)
+	AX_TICKCOLOR(ax1)	= GL_TICKCOLOR(ap)
+	AX_TICKCOLOR(ax2)	= GL_TICKCOLOR(ap)
+	AX_GRIDCOLOR(ax1)	= GL_GRIDCOLOR(ap)
+	AX_GRIDCOLOR(ax2)	= GL_GRIDCOLOR(ap)
+	AX_AXISLABELSIZE(ax1)	= GL_AXISLABELSIZE(ap)
+	AX_AXISLABELSIZE(ax2)	= GL_AXISLABELSIZE(ap)
+	AX_AXISLABELCOLOR(ax1)	= GL_AXISLABELCOLOR(ap)
+	AX_AXISLABELCOLOR(ax2)	= GL_AXISLABELCOLOR(ap)
+	AX_AXISWIDTH(ax1)	= GL_AXISWIDTH(ap)
+	AX_AXISWIDTH(ax2)	= GL_AXISWIDTH(ap)
+	AX_AXISCOLOR(ax1)	= GL_AXISCOLOR(ap)
+	AX_AXISCOLOR(ax2)	= GL_AXISCOLOR(ap)
+	AX_MINORWIDTH(ax1)	= GL_MINORWIDTH(ap)
+	AX_MINORWIDTH(ax2)	= GL_MINORWIDTH(ap)
+	AX_MAJORWIDTH(ax1)	= GL_MAJORWIDTH(ap)
+	AX_MAJORWIDTH(ax2)	= GL_MAJORWIDTH(ap)
+
+	# Compute the number of digits of precision needed for the tick labels.
+	AX_NDIGITS(ax1) = max (1, gt_ndigits (p1, p2, step))
+	AX_NDIGITS(ax2) = AX_NDIGITS(ax1)
+
+	# If both axes are to be drawn label ticks if enabled.  If only one
+	# axis is to be drawn that is the axis that must be labelled.
+
+	if (GL_DRAWAXES(ap) > 0) {
+	    AX_LABELTICKS(ax1) = GL_LABELTICKS(ap)
+	    AX_LABELTICKS(ax2) = GL_LABELTICKS(ap)
+	}
+	if (GL_DRAWAXES(ap) == 1 || GL_DRAWAXES(ap) == 3)
+	    AX_LABELTICKS(ax2) = NO
+	else if (GL_DRAWAXES(ap) == 2)
+	    AX_LABELTICKS(ax1) = NO
+
+	# The user may override the tick label format if desired.
+	if (GL_TICKFORMAT(ap) == EOS) {
+	    call sprintf (AX_TICKFORMAT(ax1), SZ_FORMAT, "%%0.%dg")
+		call pargi (AX_NDIGITS(ax1) + 1)
+	} else
+	    call strcpy (GL_TICKFORMAT(ap), AX_TICKFORMAT(ax1), SZ_FORMAT)
+	call strcpy (AX_TICKFORMAT(ax1), AX_TICKFORMAT(ax2), SZ_FORMAT)
+end
+
+
+# GLB_MINORSTEP -- Determine the step size for the minor ticks.  Adapted
+# from a routine by J. Eisenhamer (STScI) which was based on some MONGO code.
+
+real procedure glb_minorstep (x1, x2, nminor)
+
+real	x1, x2		#I interval between major ticks
+int	nminor		#I suggested number of minor ticks, or actual# if neg
+
+int	iexp
+real	amant, diff, num, range
+
+begin
+	range = abs (x2 - x1)
+	if (nminor < 0)
+	    return (range / real (-nminor + 1))
+	else {
+	    # Determine magnitude of the intervals.
+	    diff = log10 (range / nminor)
+	    iexp = int (diff)
+	    if (diff < 0)
+		iexp = iexp - 1
+	    amant = diff - real(iexp)
+
+	    # Determine an appropriate step size.
+	    if (amant < 0.15)
+		num = 1.0
+	    else if (amant < 0.50)
+		num = 2.0
+	    else if (amant < 0.85)
+		num = 5.0
+	    else
+		num = 10.0
+
+	    return (num * 10.0**iexp)
+	}
+end
diff --git a/sys/gio/glabax/glbgrid.x b/sys/gio/glabax/glbgrid.x
new file mode 100644
index 00000000..ecb24ffb
--- /dev/null
+++ b/sys/gio/glabax/glbgrid.x
@@ -0,0 +1,54 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gset.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_DRAWGRID -- Draw a grid across the plotting surface, i.e., draw
+# dotted lines between the major tick marks.
+
+procedure glb_drawgrid (gp, ax1, ax2)
+
+pointer	gp			# graphics descriptor
+pointer	ax1			# descriptor for first axis
+pointer	ax2			# descriptor for second axis
+
+int	wcs, major_tick
+real	x, y, tolerance
+real	x1, y1, x2, y2, sx, sy
+int	glb_gettick()
+errchk	glb_gettick, gseti, gsetr, gline, gctran
+
+begin
+	tolerance = TOL
+	wcs = GP_WCS(gp)
+
+	# Cache the NDC coordinates of the ends of an axis.
+	call gctran (gp, AX_START(ax1,1), AX_START(ax1,2), x1,y1, wcs, 0)
+	call gctran (gp,   AX_END(ax1,1),   AX_END(ax1,2), x2,y2, wcs, 0)
+
+	# Set polyline linetype for a dotted line.
+	call gseti (gp, G_PLTYPE, GL_DOTTED)
+	call gsetr (gp, G_PLWIDTH, 1.0)
+
+	AX_NLEFT(ax1) = -1
+	while (glb_gettick (gp, ax1, x, y, major_tick) != EOF) {
+	    if (major_tick == NO)
+		next
+
+	    # Draw grid line if we are at a major tick, provided the tick
+	    # is not at the end of the axis.
+
+	    call gctran (gp, x,y, sx,sy, wcs, 0)
+	    if (AX_HORIZONTAL(ax1) == YES) {
+		if (sx - x1 > tolerance && sx - x2 < tolerance)
+		    call gline (gp, x, AX_END(ax1,2), x, AX_END(ax2,2))
+	    } else {
+		if (sy - y1 > tolerance && sy - y2 < tolerance)
+		    call gline (gp, AX_END(ax1,1), y, AX_END(ax2,1), y)
+	    }
+	}
+
+	call gseti (gp, G_PLTYPE, GL_SOLID)
+end
diff --git a/sys/gio/glabax/glbgtick.x b/sys/gio/glabax/glbgtick.x
new file mode 100644
index 00000000..cc70fd3a
--- /dev/null
+++ b/sys/gio/glabax/glbgtick.x
@@ -0,0 +1,252 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_GETTICK -- Get the position and type of the next tick on an axis.
+# Ticks are accessed sequentially.  There are three types of tick scalings,
+# LINEAR, LOG, and ELOG.  The tick scaling need not necessarily agree with
+# the WCS scaling, hence linear tick scaling might be used on a nonlinear
+# coordinate system.  If the scaling is linear then the first tick need not
+# fall at the endpoint of the axis.  If log (or elog) scaling is in use then
+# the axis will have been rounded out to a decade and the first tick will
+# necessarily fall on the axis endpoint.  The scalings are described by the
+# following parameters:
+#
+#			(variables)
+#	nleft			number of minor ticks left to next major tick
+#	step(x,y)		displacement between minor ticks (world coords)
+#
+#			(constants)
+#	tick1(x,y)		world coords of the first tick on an axis
+#	nminor			number of minor ticks between major ticks
+#	istep(x,y)		initial step (actual step may differ)
+#	kstep(x,y)		adjustment to step at major ticks
+# 
+# KSTEP is unity if the scaling is linear.  The log scalings have a KSTEP of
+# either 10.0 or 0.1.  A negative KSTEP value is used to flag ELOG scaling.
+# ELOG, or extended range log scaling, is a log scaling which is defined for
+# X <=0 as well as x > 0.  This function is logarithmic for values less than
+# -10 or greater than 10, and linear in the range [-10:+10].  This complicates
+# tick computation because the usual 8 minor ticks per decade characteristic
+# of log scaling are not appropriate in the linear regime.  If the scaling
+# is ELOG then we ignore NMINOR and ISTEP in the linear range, changing the
+# values of these parameters temporarily to reflect 4 minor ticks with a tick
+# spacing of 2.0.
+#
+# (Note to the reader: don't feel discouraged if you don't understand this
+# (stuff, it is so complicated I don't understand it either!  Having to deal
+# (with linear, log, and elog scaling with both major and minor ticks,
+# (sometimes no minor ticks, with the axis starting at any part of the scale
+# (seems an inherently difficult problem to program compactly.  Barring
+# (programming each case separately, the best approach I could come up with was
+# (to walkthrough the code separately for each case, from all initial
+# (conditions, until it works for all cases.  If you have problems determine
+# (the initial conditions (the case) and do a similar walkthough.  Of course,
+# (if you make a change affecting one case, you may well make the code fail for
+# (a different case.
+
+int procedure glb_gettick (gp, ax, x, y, major_tick)
+
+pointer	gp			# graphics descriptor
+pointer	ax			# axis descriptor
+real	x, y			# coordinates of next tick (output)
+int	major_tick		# YES if next tick is a major tick
+
+int	i, axis, wcs, w, scaling, nminor, expon
+real	kstep, step, astep, ten, sx, sy, tolerance, pos, norm_pos
+bool	glb_eq()
+define	logscale_ 91
+
+begin
+	if (AX_DRAWTICKS(ax) == NO)
+	    return (EOF)
+
+	tolerance = TOL
+	scaling = AX_SCALING(ax)
+	nminor  = AX_NMINOR(ax)
+	kstep   = AX_KSTEP(ax)
+
+	if (AX_HORIZONTAL(ax) == YES)
+	    axis = 1
+	else
+	    axis = 2
+
+	# Count down a minor tick.  If nleft is negative then we are being
+	# called for the first time for this axis.
+
+	if (AX_NLEFT(ax) < 0) {
+
+	    # Initialize everything and return coords of the first tick.
+	    AX_KSTEP(ax) = AX_IKSTEP(ax)
+	    AX_NLEFT(ax) = AX_INLEFT(ax)
+	    do i = 1, 2 {
+		AX_POS(ax,i)  = AX_TICK1(ax,i)
+		AX_STEP(ax,i) = AX_ISTEP(ax,i)
+	    }
+
+	    step  = AX_STEP(ax,axis)
+	    astep = abs (step)
+
+	    if (AX_NLEFT(ax) == 0) {
+		# Note that there may not be any minor ticks.
+		major_tick = YES
+		AX_NLEFT(ax) = nminor
+		if (nminor > 0)
+		    if (scaling == ELOG && (astep >= .99 && astep < 2.0)) {
+			# Elog scaling in linear region.
+			AX_NLEFT(ax) = 4
+			if (step < 0)
+			    step = -2.0
+			else
+			    step =  2.0
+			AX_STEP(ax,axis) = step
+		    }
+	    } else {
+		AX_NLEFT(ax) = AX_NLEFT(ax) - 1
+		major_tick = NO
+	    }
+
+	    # Elog scaling in linear region.  KSTEP must be inverted as we
+	    # pass through the origin.  This normally occurs upon entry to the
+	    # linear region, but if we start out at +/- 10 we must set KSTEP
+	    # to its linear value during setup.
+
+	    if (scaling == ELOG && glb_eq(step,2.0))
+		AX_KSTEP(ax) = -10.0
+
+	} else {
+	    # All ticks after the first tick.
+	    do i = 1, 2
+		AX_POS(ax,i) = AX_POS(ax,i) + AX_STEP(ax,i)
+	    AX_NLEFT(ax) = AX_NLEFT(ax) - 1
+
+	    # If we are log scaling the ticks will never have more than 2
+	    # digits of precision.  Try to correct for the accumulation of
+	    # error by rounding.  When log scaling the error increases by
+	    # a factor of ten in each decade and can get quite large if
+	    # the log scale covers a large range.
+
+	    if (scaling != LINEAR) {
+		pos = AX_POS(ax,axis)
+		call fp_normr (pos, norm_pos, expon)
+		pos = nint (norm_pos * 10.0) / 10.0
+		pos = pos * (10.0 ** expon)
+		AX_POS(ax,axis) = pos
+	    }
+
+	    if (AX_NLEFT(ax) < 0) {
+		# Next tick is a major tick.  If log scaling we must reset
+		# the tick parameters for the next decade.
+
+		major_tick = YES
+		AX_NLEFT(ax) = nminor
+
+		# The following handles the special case of ELOG scaling in
+		# the linear regime when the number of minor ticks is zero.
+		# The step size in such a case is 9 to some power in the log
+		# region and +/- 10 in the linear region.
+
+		if (scaling == ELOG && nminor == 0) {
+		    pos = AX_POS(ax,axis)
+		    if (step < 0)
+			ten = -10.
+		    else
+			ten = 10.
+
+		    if (glb_eq (pos, 10.0)) {
+			if (glb_eq (step, 10.0)) {
+			    if (step < 0)
+				AX_STEP(ax,axis) = -9.
+			    else
+				AX_STEP(ax,axis) = 9.
+			    goto logscale_
+			} else
+			    step = ten
+		    } else if (glb_eq (pos, 0.0)) {
+			step = ten
+			if (pos / step < 0)
+			    AX_KSTEP(ax) = -0.1
+			else
+			    AX_KSTEP(ax) = -10.0
+		    } else
+			goto logscale_
+		    AX_STEP(ax,axis) = step
+
+		} else if (scaling != LINEAR) {
+		    # Adjust the tick step by the kstep factor, provided we
+		    # are not at the origin in ELOG scaling (the step is 1
+		    # on either side of the origin for ELOG scaling).  Reset
+		    # the step size to 1.0 if ELOG scaling and just coming out
+		    # of the linear regime.
+logscale_
+		    step = AX_STEP(ax,axis)
+		    if (scaling != ELOG || abs(AX_POS(ax,axis)) > 0.1) {
+			if (scaling == ELOG && glb_eq (step, 2.0))
+			    AX_STEP(ax,axis) = step / 2.0
+
+			do i = 1, 2
+			    AX_STEP(ax,i) = AX_STEP(ax,i) * abs (AX_KSTEP(ax))
+		    }
+
+		    # Adjust the step size to 2.0 if ELOG scaling and in the
+		    # linear regime (initial step size of 1).
+
+		    step = AX_STEP(ax,axis)
+		    if (scaling == ELOG && glb_eq(step,1.0)) {
+			if (step < 0)
+			    step = -2.0
+			else
+			    step =  2.0
+			AX_STEP(ax,axis) = step
+		    }
+
+		    # If elog scaling and we have just entered the linear
+		    # regime, adjust the number of ticks and the KSTEP factor.
+
+		    if (scaling == ELOG && glb_eq(step,2.0)) {
+			# Elog scaling in linear region.  KSTEP must be
+			# inverted as we pass through the origin.
+
+			if (abs(AX_POS(ax,axis)) > 0.1)
+			    AX_KSTEP(ax) = -10.0
+
+			if (nminor > 0)
+			    AX_NLEFT(ax) = 4
+		    }
+		}
+	    } else
+		major_tick = NO
+	}
+
+	x = AX_POS(ax,1)
+	y = AX_POS(ax,2)
+
+	# Return EOF if tick falls beyond end of axis.  The comparison is made
+	# in NDC coords to avoid having to check if the WCS is increasing or
+	# decreasing and to avoid the problems of comparing unnormalized
+	# floating point numbers.
+
+	wcs = GP_WCS(gp)
+	w = GP_WCSPTR(gp,wcs)
+
+	call gctran (gp, x,y, sx,sy, wcs, 0)
+	if (sx - WCS_SX2(w) > tolerance || sy - WCS_SY2(w) > tolerance)
+	    return (EOF)
+	else
+	    return (OK)
+end
+
+
+# GLB_EQ -- Compare two (near normalized) floating point numbers for
+# equality, using the absolute value of the first argument.
+
+bool procedure glb_eq (a, b)
+
+real	a			# compare absolute value of this number
+real	b			# to this positive number
+
+begin
+	return (abs (abs(a) - b) < 0.1)
+end
diff --git a/sys/gio/glabax/glblabel.x b/sys/gio/glabax/glblabel.x
new file mode 100644
index 00000000..ecf57c94
--- /dev/null
+++ b/sys/gio/glabax/glblabel.x
@@ -0,0 +1,84 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gset.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_LABEL_AXIS -- Label an axis.  If both axes were drawn only the first is
+# labelled, otherwise the label is placed on withever axis was drawn.  This is
+# done by drawing the axis labels just outside the tick mark labels, wherever
+# those happened to be.  The axis label offset is in the same direction as the
+# tick label offset and is centered on each axis.  The distance from the axis
+# is a function of the size of the tick labels.
+
+procedure glb_label_axis (gp, ax, xlabel, ylabel)
+
+pointer	gp			# graphics descriptor
+pointer	ax			# axis descriptor
+char	xlabel[ARB]		# X axis label
+char	ylabel[ARB]		# Y axis label
+
+int	wcs
+real	x1, x2, y1, y2, x, y, dx, dy
+real	char_height, char_width
+int	strlen()
+real	ggetr()
+
+begin
+	wcs = GP_WCS(gp)
+
+	# Get character height and width in NDC coords.
+	char_height = ggetr (gp, "ch")
+	char_width  = ggetr (gp, "cw")
+
+	if (char_height < EPSILON)
+	    char_height = DEF_CHARHEIGHT
+	if (char_width < EPSILON)
+	    char_width = DEF_CHARWIDTH
+
+	# Compute axis center in NDC coords.
+	call gctran (gp, AX_START(ax,1), AX_START(ax,2), x1,y1, wcs, 0)
+	call gctran (gp,   AX_END(ax,1),   AX_END(ax,2), x2,y2, wcs, 0)
+	x = (x1 + x2) / 2.0
+	y = (y1 + y2) / 2.0
+
+	# Set relative text size and get device character size for a text
+	# size of 1.0.  Set WCS to NDC coords since the offset to the
+	# tick label is in NDC coordinates.
+
+	call gsetr (gp, G_TXSIZE, AX_AXISLABELSIZE(ax))
+	call gseti (gp, G_WCS, 0)
+
+	# Draw the axis label.
+
+	if (AX_HORIZONTAL(ax) == YES) {
+	    # Axis is horizontal.  Tick label vector tells us whether to
+	    # draw axis label above or below axis.
+
+	    if (strlen (xlabel) > 0) {
+		dy = 2.0 * AX_TICKLABELSIZE(ax) * char_height +
+		     0.5 * AX_AXISLABELSIZE(ax) * char_height
+		if (AX_TICKLABELOFFSET(ax,2) < 0)
+		    dy = -dy
+		call gtext (gp, x, y + dy, xlabel, "hj=c;vj=c")
+	    }
+	} else {
+	    # Axis is vertical.  Always put label fixed distance from axis
+	    # regardless of size of tick labels (for consistency and to
+	    # avoid clipping at the device screen boundary).  Label runs
+	    # bottom to top in a vertical field with char up pointing to
+	    # the left.
+
+	    if (strlen (ylabel) > 0) {
+		dx = (Y_LABELOFFSET * char_width * AX_TICKLABELSIZE(ax)) + 
+		    0.5 * AX_AXISLABELSIZE(ax) * char_height
+		    
+		if (AX_TICKLABELOFFSET(ax,1) < 0)
+		    dx = -dx
+		call gtext (gp, x + dx, y, ylabel, "up=180;hj=c;vj=c")
+	    }
+	}
+
+	call gseti (gp, G_WCS, wcs)
+end
diff --git a/sys/gio/glabax/glbloglab.x b/sys/gio/glabax/glbloglab.x
new file mode 100644
index 00000000..6e7ec1cc
--- /dev/null
+++ b/sys/gio/glabax/glbloglab.x
@@ -0,0 +1,139 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gset.h>
+include	<gio.h>
+include	"glabax.h"
+
+define	SZ_MANTISSA	3		# "10" or "-10"
+define	SZ_EXPONENT	4		# largest is "-999"
+
+
+# GLB_LOGLAB -- Draw a tick label in log units at the given position.
+# A log tick is a power of ten, e.g.  10^2, where the ^ signifies that
+# the 2 is to be drawn one half character height higher than the 10.
+
+procedure glb_loglab (gp, sx, sy, val, fmt, scaling)
+
+pointer	gp			# graphics descriptor
+real	sx, sy			# NDC coords of label
+real	val			# value to be encoded (not the log of)
+char	fmt[ARB]		# tick label gtext format (justification)
+int	scaling			# type of scaling on axis
+
+bool	zero
+char	mantissa[SZ_MANTISSA]
+char	exponent[SZ_EXPONENT]
+int	len_mantissa, len_exponent, ip, hj, vj
+real	logval, char_height, char_width, left, xpos, ypos, txsize
+
+bool	fp_equalr()
+real	elogr(), gstatr(), ggetr()
+int	strlen(), strmatch(), itoc()
+
+begin
+	# Compute the log value to be encoded.
+	if (scaling == LOG)
+	    logval = log10 (val)
+	else {
+	    logval = elogr (val)
+	    zero = fp_equalr (logval, 0.0)
+	}
+
+	txsize = gstatr (gp, G_TXSIZE)
+
+	# Get char height and width in NDC coords.
+	char_height = ggetr (gp, "ch")
+	if (char_height < EPSILON)
+	    char_height = DEF_CHARHEIGHT
+	char_height = char_height * txsize
+
+	char_width = ggetr (gp, "cw")
+	if (char_width < EPSILON)
+	    char_width = DEF_CHARWIDTH
+	char_width = char_width * txsize
+
+	# Encode the mantissa and exponent strings.
+	if (zero) {
+	    call strcpy ("0", mantissa, SZ_MANTISSA)
+	} else if (logval < 0 && scaling == ELOG) {
+	    call strcpy ("-10", mantissa, SZ_MANTISSA)
+	    logval = abs (logval)
+	} else
+	    call strcpy ("10", mantissa, SZ_MANTISSA)
+
+	len_mantissa = strlen (mantissa)
+	if (zero)
+	    len_exponent = 0
+	else
+	    len_exponent = itoc (nint(logval), exponent, SZ_EXPONENT)
+
+	# Determine type of horizontal justification required.
+	ip = strmatch (fmt, "hj=")
+	if (ip <= 0)
+	    hj = 'c'
+	else
+	    hj = fmt[ip]
+
+	# Determine type of vertical justification required.
+	ip = strmatch (fmt, "vj=")
+	if (ip <= 0)
+	    vj = 'c'
+	else
+	    vj = fmt[ip]
+
+	# On devices with adjustable character sizes the most pleasing results
+	# are obtained if the digits "10" are nicely aligned on the vertical
+	# axis, regardless of the actual number of characters in the exponent
+	# string, minus signs etc (this type of alignment is more natural
+	# because the exponent is printed at half size).  Hence if we are on
+	# a vertical axis (hj != c) fix the number of characters in the two
+	# strings so that the alignment comes out the same regardless of the
+	# actual number of chars in either field.  The length of the exponent
+	# field is not completely fixed, rather we allow a little more space
+	# if the exponent is large.  For small exponents len_exponent=1.
+
+	if (hj != 'c') {
+	    len_mantissa = 2
+	    len_exponent = (len_exponent + 1) / 2
+	}
+
+	# Compute XPOS, the NDC X coord of the point halfway between the
+	# last char of the mantissa and the first char of the exponent.
+
+	switch (hj) {
+	case 'l':
+	    left = sx
+	case 'r':
+	    left = sx - (len_mantissa + len_exponent) * char_width
+	default:
+	    left = sx - ((len_mantissa + len_exponent) * char_width) / 2.0
+	}
+
+	xpos = left + len_mantissa * char_width
+
+	# Compute YPOS, the NDC Y coord of the center of a mantissa character
+	# and of the bottom of an exponent character.  Using the same coordinate
+	# to address both positions makes the label come out the same regardless
+	# of the plot magnification, even on a device where the character size
+	# is fixed by the hardware.
+
+	switch (vj) {
+	case 'b':
+	    ypos = sy + char_height / 2.0
+	case 't':
+	    ypos = sy - char_height / 2.0
+	default:
+	    ypos = sy
+	}
+
+	# Draw the mantissa.
+	call gtext (gp, xpos, ypos, mantissa, "hj=r,vj=c")
+
+	# Draw the exponent if there is one.
+	if (!zero) {
+	    call gsetr (gp, G_TXSIZE, txsize / 2.0)
+	    call gtext (gp, xpos, ypos, exponent, "hj=l;vj=b")
+	    call gsetr (gp, G_TXSIZE, txsize)
+	}
+end
diff --git a/sys/gio/glabax/glbsetax.x b/sys/gio/glabax/glbsetax.x
new file mode 100644
index 00000000..f0c9aa29
--- /dev/null
+++ b/sys/gio/glabax/glbsetax.x
@@ -0,0 +1,130 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_SET_AXES -- Set all axis descriptor parameters not pertaining to the
+# ticks.  The WCS has already been fixed by the time we get here.
+
+procedure glb_set_axes (gp, ap, ax1, ax2, angle)
+
+pointer	gp			# graphics descriptor
+pointer	ap			# axis parameters (from graphics descriptor)
+pointer	ax1, ax2		# axis descriptors (output)
+int	angle			# axis orientation, 0 or 90 degrees
+
+pointer	w
+int	axis
+real	p1, p2
+real	x1, x2, y1, y2
+real	glb_ticklen()
+
+begin
+	w = GP_WCSPTR (gp, GP_WCS(gp))
+
+	# If the window was rounded in Y in the second call to find_ticks,
+	# then the Y positions of the first ticks set in the first call will
+	# be in error and must be corrected.  If the user has elected to set
+	# the axis position explicitly, however, then we must leave it alone.
+
+	if (angle == 0 && GL_SETAXISPOS(GP_XAP(gp)) == NO) {
+	    AX_TICK1(ax1,2) = WCS_WY1(w)
+	    AX_TICK1(ax2,2) = WCS_WY2(w)
+	}
+
+	# Set the tick lengths.  This is done here rather than in findticks
+	# due to rounding, as noted above.  The tick offsets in world
+	# coordinates.  The GL values are given in NDC coordinates.
+
+	if (angle == 0) {
+	    axis = 2
+	    AX_HORIZONTAL(ax1) = YES
+	    AX_HORIZONTAL(ax2) = YES
+	} else {
+	    axis = 1
+	    AX_HORIZONTAL(ax1) = NO
+	    AX_HORIZONTAL(ax2) = NO
+	}
+
+	AX_MAJORTICK(ax1,axis) = glb_ticklen (gp, ax1,  GL_MAJORLENGTH(ap))
+	AX_MINORTICK(ax1,axis) = glb_ticklen (gp, ax1,  GL_MINORLENGTH(ap))
+	AX_MAJORTICK(ax2,axis) = glb_ticklen (gp, ax2, -GL_MAJORLENGTH(ap))
+	AX_MINORTICK(ax2,axis) = glb_ticklen (gp, ax2, -GL_MINORLENGTH(ap))
+
+	# Select none, either, or both axes to be drawn.  If only the second
+	# axis is drawn then that is the side we must draw the tick and axis
+	# labels on.
+
+	switch (GL_DRAWAXES(ap)) {
+	case 0:
+	    AX_DRAWME(ax1) = NO
+	    AX_DRAWME(ax2) = NO
+	    return
+	case 1:
+	    AX_DRAWME(ax1) = YES
+	    AX_DRAWME(ax2) = NO
+	case 2:
+	    AX_DRAWME(ax1) = NO
+	    AX_DRAWME(ax2) = YES
+	default:
+	    AX_DRAWME(ax1) = YES
+	    AX_DRAWME(ax2) = YES
+	}
+
+	# Determine the endpoints of the axis.  These default to the corners of
+	# the viewport (in world coordinates), but the positions may be
+	# overriden by the user if desired.
+
+	# First get the positions of the two axes.
+	if (GL_SETAXISPOS(ap) == YES) {
+	    p1 = GL_AXISPOS1(ap)
+	    p2 = GL_AXISPOS2(ap)
+	} else if (angle == 0) {
+	    p1 = WCS_WY1(w)
+	    p2 = WCS_WY2(w)
+	} else {
+	    p1 = WCS_WX1(w)
+	    p2 = WCS_WX2(w)
+	}
+
+	# Convert these positions into the world coordinates of the endpoints.
+	if (angle == 0) {
+	    x1 = WCS_WX1(w)
+	    x2 = WCS_WX2(w)
+	    y1 = p1
+	    y2 = p2
+	} else {
+	    x1 = p1
+	    x2 = p2
+	    y1 = WCS_WY1(w)
+	    y2 = WCS_WY2(w)
+	}
+
+	if (angle == 0) {
+	    # Set the left and right endpoints of the axes.
+
+	    AX_START(ax1,1) = x1
+	    AX_START(ax1,2) = y1
+	    AX_END(ax1,1)   = x2
+	    AX_END(ax1,2)   = y1
+
+	    AX_START(ax2,1) = x1
+	    AX_START(ax2,2) = y2
+	    AX_END(ax2,1)   = x2
+	    AX_END(ax2,2)   = y2
+
+	} else {
+	    # Set the lower and upper endpoints of the axes.
+
+	    AX_START(ax1,1) = x1
+	    AX_START(ax1,2) = y1
+	    AX_END(ax1,1)   = x1
+	    AX_END(ax1,2)   = y2
+
+	    AX_START(ax2,1) = x2
+	    AX_START(ax2,2) = y1
+	    AX_END(ax2,1)   = x2
+	    AX_END(ax2,2)   = y2
+	}
+end
diff --git a/sys/gio/glabax/glbsetup.x b/sys/gio/glabax/glbsetup.x
new file mode 100644
index 00000000..a609d2ad
--- /dev/null
+++ b/sys/gio/glabax/glbsetup.x
@@ -0,0 +1,51 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<gio.h>
+
+# GLB_SETUP -- Set up the axis drawing and labelling parameters.  These are
+# the coordinate transformations, i.e., log scaling, window and viewport
+# coordinates, plus the parameters which pertain only to axis drawing and
+# labelling.  The order in which the subprocedures are called is significant.
+
+procedure glb_setup (gp, axes, ntitlelines, xlabel, ylabel)
+
+pointer	gp			# graphics descriptor
+pointer	axes[4]			# array of pointers to axis descriptors
+int	ntitlelines		# number of lines in title block
+char	xlabel[ARB]		# x axis label
+char	ylabel[ARB]		# y axis label
+
+pointer	w
+bool	fp_nondegenr()
+
+begin
+	w = GP_WCSPTR (gp, GP_WCS(gp))
+
+	# Verify that there is sufficient range in the wcs X and Y.
+	if (fp_nondegenr (WCS_WX1(w), WCS_WX2(w)))
+	    GP_WCSSTATE(gp) = MODIFIED
+	if (fp_nondegenr (WCS_WY1(w), WCS_WY2(w)))
+	    GP_WCSSTATE(gp) = MODIFIED
+
+	# If log scaling is in effect on either axis, verify that log scaling
+	# is sensible and if so select either LOG or ELOG scaling.
+
+	call glb_verify_log_scaling (gp)
+
+	# Set the viewport if not already set.
+	call glb_set_viewport (gp, ntitlelines, xlabel, ylabel)
+
+	# Find the best positions for the tick marks, and if rounding is
+	# enabled, extend the WCS outward to the next tick mark on either
+	# end.
+
+	call glb_find_ticks (gp, GP_XAP(gp), axes[1], axes[4],  0)
+	call glb_find_ticks (gp, GP_YAP(gp), axes[3], axes[2], 90)
+
+	# Set the remaining parameters in the axis drawing descriptors.
+	# Must not be called until the window and viewport coordinates are
+	# fixed.
+
+	call glb_set_axes (gp, GP_XAP(gp), axes[1], axes[4],  0)
+	call glb_set_axes (gp, GP_YAP(gp), axes[3], axes[2], 90)
+end
diff --git a/sys/gio/glabax/glbsview.x b/sys/gio/glabax/glbsview.x
new file mode 100644
index 00000000..1b099b1a
--- /dev/null
+++ b/sys/gio/glabax/glbsview.x
@@ -0,0 +1,117 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_SET_VIEWPORT -- If the viewport has not yet been set, i.e., if the
+# viewport is still [0:1,0:1], compute the size of the largest viewport which
+# leaves sufficient room around the border for the axis labels and plot title.
+# If a nonzero aspect ratio is specified make the viewport have that aspect
+# ratio.
+
+procedure glb_set_viewport (gp, ntitlelines, xlabel, ylabel)
+
+pointer	gp			# graphics descriptor
+int	ntitlelines		# number of lines to reserve for title block
+char	xlabel[ARB]		# x axis label
+char	ylabel[ARB]		# y axis label
+
+pointer	w, xap, yap
+bool	draw_title, draw_xlabel, draw_ylabel, draw_xticks, draw_yticks
+real	char_height, char_width
+real	aspect, cur_aspect, dev_aspect, dx, dy
+real	xwidth, ywidth, yreserve
+real	ggetr()
+
+begin
+	w = GP_WCSPTR (gp, GP_WCS(gp))
+	xap = GP_XAP(gp)
+	yap = GP_YAP(gp)
+
+	if ((WCS_SX1(w) > EPSILON) || (abs(1.0 - WCS_SX2(w)) > EPSILON) ||
+	    (WCS_SY1(w) > EPSILON) || (abs(1.0 - WCS_SY2(w)) > EPSILON))
+	    return
+
+	draw_title  = (ntitlelines > 0 && GP_DRAWTITLE(gp) == YES)
+	draw_xticks = (GL_DRAWAXES(xap) > 0 && GL_LABELTICKS(xap) == YES)
+	draw_xlabel =
+	    (draw_xticks && xlabel[1] != EOS && GL_LABELAXIS(xap) == YES)
+	draw_yticks = (GL_DRAWAXES(yap) > 0 && GL_LABELTICKS(yap) == YES)
+	draw_ylabel =
+	    (draw_yticks && ylabel[1] != EOS && GL_LABELAXIS(yap) == YES)
+
+	char_width  = ggetr (gp, "cw")
+	char_height = ggetr (gp, "ch")
+
+	if (char_width < EPSILON)
+	    char_width = DEF_CHARWIDTH
+	if (char_height < EPSILON)
+	    char_height = DEF_CHARHEIGHT
+
+	# X axis.
+	if (draw_yticks && draw_ylabel)
+	    xwidth = max (4, LEFT_BORDER + 2)
+	else if (draw_yticks)
+	    xwidth = max (4, LEFT_BORDER)
+	else
+	    xwidth = 0
+	xwidth = xwidth * char_width * GL_TICKLABELSIZE(xap)
+	    
+	# Y axis.
+	if (draw_xticks && draw_xlabel)
+	    ywidth = BOTTOM_BORDER
+	else if (draw_xticks)
+	    ywidth = max (2, (BOTTOM_BORDER - 2))
+	else
+	    ywidth = 0
+	ywidth = ywidth * char_height * GL_TICKLABELSIZE(yap)
+
+	# Compute amount of extra space to allow for the title block, which
+	# may contain more than one line.
+
+	if (!draw_title && !draw_xticks && !draw_yticks)
+	    yreserve = 0
+	else if (!draw_title && GP_ASPECT(gp) > 0.9)
+	    yreserve = 0
+	else {
+	    yreserve = min (MAX_SZTITLEBLOCK,
+		max (MIN_NTITLELINES, ntitlelines + 1) *
+		char_height * GP_TITLESIZE(gp))
+	}
+
+	# Set the viewport.  The viewport is the largest area yielding the
+	# desired borders.  The viewport is centered in X and positioned just
+	# below the title block in Y.
+
+	WCS_SX1(w) = xwidth
+	WCS_SX2(w) = 1.0 - xwidth
+	WCS_SY1(w) = ywidth
+	WCS_SY2(w) = 1.0 - yreserve
+
+	# Adjust the viewport to achieve the specified aspect ratio, if a
+	# nonzero aspect ratio was given.
+
+	dev_aspect = GP_DEVASPECT(gp)		# device aspect ratio
+	aspect     = GP_ASPECT(gp)		# user desired aspect ratio
+
+	if (aspect > EPSILON) {
+	    dx = WCS_SX2(w) - WCS_SX1(w)
+	    dy = WCS_SY2(w) - WCS_SY1(w)
+	    cur_aspect = dy / dx * dev_aspect
+
+	    if (cur_aspect > aspect) {
+		# Viewport is taller than desired.
+		dy = aspect / dev_aspect * dx
+		WCS_SY1(w) = (1.0 - dy) / 2.0
+		WCS_SY2(w) = 1.0 - WCS_SY1(w)
+	    } else {
+		# Viewport is not as wide as desired.
+		dx = dev_aspect * dy / aspect
+		WCS_SX1(w) = (1.0 - dx) / 2.0
+		WCS_SX2(w) = 1.0 - WCS_SX1(w)
+	    }
+	}
+
+	GP_WCSSTATE(gp) = MODIFIED
+end
diff --git a/sys/gio/glabax/glbticlen.x b/sys/gio/glabax/glbticlen.x
new file mode 100644
index 00000000..de557757
--- /dev/null
+++ b/sys/gio/glabax/glbticlen.x
@@ -0,0 +1,42 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_TICKLEN -- Compute the length of a tick in world coordinates.  All tick
+# drawing is performed in world coordinates since the ticks show the world
+# coordinate system.  The position of a tick must be computed in world coords
+# when the axis is drawn to reflect log scaling (or any other nonlinear
+# scaling.  Less obviously, the tick offset should be given in world coords
+# so that when the tick is drawn by a GRDRAW the tick will follow a line of
+# constant X or Y in world coordinates, and this line will not necessarily be
+# a line of constant X or Y in NDC coordinates.
+
+real procedure glb_ticklen (gp, ax, ndc_length)
+
+pointer	gp			# graphics descriptor
+pointer	ax			# axis descriptor
+real	ndc_length		# length of tick in NDC units
+
+int	wcs
+real	x, y, wx, wy
+
+begin
+	wcs = GP_WCS(gp)
+	call gctran (gp, AX_TICK1(ax,1), AX_TICK1(ax,2), x, y, wcs, 0)
+
+	if (AX_HORIZONTAL(ax) == YES)
+	    y = y + ndc_length
+	else
+	    x = x + ndc_length
+
+	call gctran (gp, x, y, wx, wy, 0, wcs)
+	if (AX_HORIZONTAL(ax) == YES) {
+	    call pargr (wy - AX_TICK1(ax,2))
+	    return (wy - AX_TICK1(ax,2))
+	} else {
+	    call pargr (wx - AX_TICK1(ax,1))
+	    return (wx - AX_TICK1(ax,1))
+	}
+end
diff --git a/sys/gio/glabax/glbtitle.x b/sys/gio/glabax/glbtitle.x
new file mode 100644
index 00000000..d8c43c67
--- /dev/null
+++ b/sys/gio/glabax/glbtitle.x
@@ -0,0 +1,76 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gset.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_PLOT_TITLE -- Draw plot title block.  The block may contain several lines.
+# Lines are plotted with center, left, or right justification, immediately
+# above the top viewport boundary (not immediately above the drawn axis,
+# which need not be at the viewport boundary).
+
+procedure glb_plot_title (gp, title, ntitlelines)
+
+pointer	gp			# graphics descriptor
+char	title[ARB]		# title block
+int	ntitlelines		# number of lines in title block
+
+int	lineno, ip, wcs
+real	char_height, x, y, dy
+pointer sp, op, lbuf, format, w
+real	ggetr()
+
+begin
+	if (title[1] == EOS || ntitlelines < 1)
+	    return
+
+	call smark (sp)
+	call salloc (lbuf, SZ_LINE, TY_CHAR)
+	call salloc (format, SZ_FORMAT, TY_CHAR)
+
+	char_height = ggetr (gp, "ch")
+	if (char_height < EPSILON)
+	    char_height = DEF_CHARHEIGHT * GP_TITLESIZE(gp)
+
+	wcs = GP_WCS(gp)
+	w = GP_WCSPTR (gp, wcs)
+	y = min (1.0 - char_height,
+	    WCS_SY2(w) + (ntitlelines - 1 + 0.5) * char_height)
+
+	call sprintf (Memc[format], SZ_FORMAT, "hj=%c,vj=b")
+	    switch (GP_TITLEJUST(gp)) {
+	    case GT_LEFT:
+		call pargi ('l')
+		x = WCS_SX1(w)
+	    case GT_RIGHT:
+		call pargi ('r')
+		x = WCS_SX2(w)
+	    default:
+		call pargi ('c')
+		x = (WCS_SX1(w) + WCS_SX2(w)) / 2.0
+	    }
+
+	call gsetr (gp, G_TXSIZE, GP_TITLESIZE(gp))
+	call gseti (gp, G_WCS, 0)
+	lineno = 1
+	op = lbuf
+
+	for (ip=1;  title[ip] != EOS;  ip=ip+1)
+	    if (title[ip] == '\n' || (title[ip+1] == EOS && op > lbuf)) {
+		if (title[ip] != '\n') {
+		    Memc[op] = title[ip]
+		    op = op + 1
+		}
+		Memc[op] = EOS
+		dy = (lineno - 1) * char_height
+		call gtext (gp, x, y - dy, Memc[lbuf], Memc[format])
+		lineno = lineno + 1
+		op = lbuf
+	    } else {
+		Memc[op] = title[ip]
+		op = op + 1
+	    }
+	
+	call sfree (sp)
+end
diff --git a/sys/gio/glabax/glbverify.x b/sys/gio/glabax/glbverify.x
new file mode 100644
index 00000000..6666b06a
--- /dev/null
+++ b/sys/gio/glabax/glbverify.x
@@ -0,0 +1,36 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<gset.h>
+include	<gio.h>
+include	"glabax.h"
+
+# GLB_VERIFY_LOG_SCALING -- Verify that log scaling makes sense, i.e., that
+# the range covered by an axis compared to its distance from the origin is
+# large enough to permit log scaling.  If log scaling is reasonable check if
+# the window goes negative, and switch to ELOG scaling if such is the case.
+
+procedure glb_verify_log_scaling (gp)
+
+pointer	gp			# graphics descriptor
+pointer	w
+
+begin
+	w  = GP_WCSPTR (gp, GP_WCS(gp))
+
+	# Force ELOG scaling if any data <= 0.
+
+	if (WCS_XTRAN(w) != LINEAR)
+	    if (WCS_WX1(w) <= 0 || WCS_WX2(w) <= 0)
+		WCS_XTRAN(w) = ELOG
+
+	if (WCS_YTRAN(w) != LINEAR)
+	    if (WCS_WY1(w) <= 0 || WCS_WY2(w) <= 0)
+		WCS_YTRAN(w) = ELOG
+
+	# Set the WCS state to modified even if it wasn't.  This is safe
+	# and in any case the WCS is changed in the main glabax routine
+	# shortly after we are called.
+
+	GP_WCSSTATE(gp) = MODIFIED
+end
diff --git a/sys/gio/glabax/mkpkg b/sys/gio/glabax/mkpkg
new file mode 100644
index 00000000..c8990e1a
--- /dev/null
+++ b/sys/gio/glabax/mkpkg
@@ -0,0 +1,22 @@
+# Make the GLABAX axis drawing and labelling package.
+
+$checkout libex.a lib$
+$update   libex.a
+$checkin  libex.a lib$
+$exit
+
+libex.a:
+	glabax.x	glabax.h <gio.h> <gset.h> <mach.h>
+	glbencode.x	glabax.h <mach.h>
+	glbfind.x	glabax.h <gio.h> <mach.h>
+	glbgrid.x	glabax.h <gio.h> <gset.h> <mach.h>
+	glbgtick.x	glabax.h <gio.h> <mach.h>
+	glblabel.x	glabax.h <gio.h> <gset.h> <mach.h>
+	glbloglab.x	glabax.h <gio.h> <gset.h> <mach.h>
+	glbsetax.x	glabax.h <gio.h> <mach.h>
+	glbsetup.x	<gio.h>
+	glbsview.x	glabax.h <gio.h> <mach.h>
+	glbticlen.x	glabax.h <gio.h> <mach.h>
+	glbtitle.x	glabax.h <gio.h> <gset.h> <mach.h>
+	glbverify.x	glabax.h <gio.h> <gset.h> <mach.h>
+	;