Initial commit

1 files changed, 940 insertions, 0 deletions

diff --git a/noao/digiphot/ptools/pexamine/ptimplot.x b/noao/digiphot/ptools/pexamine/ptimplot.x
new file mode 100644
index 00000000..8e586b97
--- /dev/null
+++ b/noao/digiphot/ptools/pexamine/ptimplot.x
@@ -0,0 +1,940 @@
+include <imhdr.h>
+include <gset.h>
+include <math.h>
+include	<mach.h>
+
+
+# PT_RPLOT -- Plot the radial profile.
+
+procedure pt_rplot (gd, im, wx, wy)
+
+pointer	gd		# pointer to the graphics stream
+pointer	im		# pointer to the input image
+real	wx, wy		# radial profile coordinates
+
+int	marker_type, lenbuf, npix
+pointer	sp, radius, intensity, marker, longtitle
+real	szmarker, rin, rout, x1, x2, y1, y2, dmin, dmax
+bool	clgetb()
+int	clgeti(), strlen(), pt_radpix()
+real	clgetr()
+
+begin
+	# Check for undefined graphics stream.
+	if (gd == NULL) {
+	    call printf ("The graphics device is undefined\n")
+	    return
+	}
+
+	# Check for undefined input image.
+	if (im == NULL) {
+	    call printf ("The graphics device is undefined\n")
+	    return
+	}
+
+	# Check for an undefined center.
+	if (IS_INDEFR(wx) || IS_INDEFR(wy)) {
+	    call printf ("The radial profile plot center is undefined.\n")
+	    return
+	}
+
+	# Get the inner and outer radii.
+	rin = clgetr ("radplot.rinner")
+	rout = clgetr ("radplot.router")
+	if (rout <= rin) {
+	    call printf (
+	   "The outer radius %g is <= the inner radius %g\n")
+		call pargr (rin)
+		call pargr (rout)
+	    return
+	}
+
+	# Allocate working space.
+	call smark (sp)
+
+	# Get the data.
+	lenbuf = PI * (rin + rout + 1.0) * (rout - rin + 0.5)
+	call salloc (radius, lenbuf, TY_REAL)
+	call salloc (intensity, lenbuf, TY_REAL)
+	npix = pt_radpix (im, wx, wy, rin, rout, Memr[radius], Memr[intensity])
+	if (npix <= 0) {
+	    call printf ("The object at %g %g is off the image\n")
+		call pargr (wx)
+		call pargr (wy)
+	    call sfree (sp)
+	    return
+	}
+
+	# Clear the plotting strucuture.
+	call gclear (gd)
+
+	# Fetch the window and viewport parameters.
+	x1 = clgetr ("radplot.x1")
+	x2 = clgetr ("radplot.x2")
+	y1 = clgetr ("radplot.y1")
+	y2 = clgetr ("radplot.y2")
+	if (IS_INDEFR(x1) || IS_INDEFR(x2)) {
+	    call pt_alimr (Memr[radius], npix, dmin, dmax)
+	    if (IS_INDEFR(x1))
+		x1 = dmin
+	    if (IS_INDEFR(x2))
+		x2 = dmax
+	}
+	if (IS_INDEFR(y1) || IS_INDEFR(y2)) {
+	    call pt_alimr (Memr[intensity], npix, dmin, dmax)
+	    if (IS_INDEFR(y1))
+		y1 = dmin
+	    if (IS_INDEFR(y2))
+		y2 = dmax
+	}
+
+	# Set the scale of the axes.
+	call gswind (gd, x1, x2, y1, y2)
+	if (clgetb ("radplot.logx"))
+	    call gseti (gd, G_XTRAN, GW_LOG)
+	else
+	    call gseti (gd, G_XTRAN, GW_LINEAR)
+	if (clgetb ("radplot.logy"))
+	    call gseti (gd, G_YTRAN, GW_LOG)
+	else
+	    call gseti (gd, G_YTRAN, GW_LINEAR)
+
+	# Get the x and y axes parameters.
+	if (! clgetb ("radplot.fill"))
+	    call gseti (gd, G_ASPECT, 1)
+	if (clgetb ("radplot.round"))
+	    call gseti (gd, G_ROUND, YES)
+
+	# Get the axis drawing parameters. 
+	if (clgetb ("radplot.box")) {
+
+	    # Get number of major and minor tick marks.
+	    call gseti (gd, G_XNMAJOR, clgeti ("radplot.majrx"))
+	    call gseti (gd, G_XNMINOR, clgeti ("radplot.minrx"))
+	    call gseti (gd, G_YNMAJOR, clgeti ("radplot.majry"))
+	    call gseti (gd, G_YNMINOR, clgeti ("radplot.minry"))
+
+	    # Label tick marks on axes.
+	    if (clgetb ("radplot.ticklabels"))
+	        call gseti (gd, G_LABELTICKS, YES)
+	    else
+	        call gseti (gd, G_LABELTICKS, NO)
+
+	    # Draw grid.
+	    if (clgetb ("radplot.grid"))
+	        call gseti (gd, G_DRAWGRID, YES)
+	    else
+	        call gseti (gd, G_DRAWGRID, NO)
+
+	    # Optionally draw a box around the plot.
+	    call salloc (longtitle, 2 * SZ_LINE, TY_CHAR)
+	    if (clgetb ("radplot.banner")) {
+		call sysid (Memc[longtitle], 2 * SZ_LINE)
+	        call sprintf (Memc[longtitle+strlen(Memc[longtitle])],
+		    2 * SZ_LINE, "\n%s: xc=%g  yc=%g  rinner=%g  router=%g\n%s")
+		    call pargstr (IM_HDRFILE(im))
+		    call pargr (wx)
+		    call pargr (wy)
+		    call pargr (rin)
+		    call pargr (rout)
+		    call pargstr ("Radial Profile Plot")
+	    } else {
+	        call sprintf (Memc[longtitle], 2 * SZ_LINE,
+	            "%s: xc=%g  yc=%g  rinner=%g  router=%g\n%s")
+		    call pargstr (IM_HDRFILE(im))
+		    call pargr (wx)
+		    call pargr (wy)
+		    call pargr (rin)
+		    call pargr (rout)
+		    call pargstr ("Radial Profile Plot")
+	    }
+	    call glabax (gd, Memc[longtitle], "Radial distance (pixels)",
+	        "Intensity (counts)")
+	}
+
+	# Get the marker type, the size of the marker and the linewidth.
+	call salloc (marker, SZ_FNAME, TY_CHAR)
+	call clgstr ("radplot.marker", Memc[marker], SZ_FNAME)
+	call pt_marker (Memc[marker], SZ_FNAME, marker_type)
+	if (marker_type != GM_POINT)
+	    szmarker = clgetr ("radplot.szmarker")
+	else
+	    szmarker = 0.0
+	call gsetr (gd, G_PLWIDTH, 2.0)
+
+	# Draw the points in using the deletions array.
+	call gpmark (gd, Memr[radius], Memr[intensity], npix, marker_type,
+	    szmarker, szmarker)
+	call gflush (gd)
+
+	call sfree (sp)
+end
+
+
+define	CSIZE		24
+
+# PT_SPLOT -- Draw a perspective view of a surface.  The altitude
+# and azimuth of the viewing angle are variable.
+ 
+procedure pt_splot (gd, im, x, y)
+ 
+pointer	gd		# pointer to the graphics stream
+pointer	im		# pointer to the image descriptor
+real	x, y		# the object center
+
+int	nx, ny, x1, x2, y1, y2, npts, wkid
+pointer	data, sp, sdata, work, longtitle
+real	floor, ceiling, angv, angh
+bool	clgetb()
+int 	clgeti(), strlen()
+pointer	pt_gdata()
+real	clgetr()
+
+int	first
+real	vpx1, vpx2, vpy1, vpy2
+common  /frstfg/ first
+common  /noaovp/ vpx1, vpx2, vpy1, vpy2
+ 
+begin
+	# Check for undefined graphics stream.
+	if (gd == NULL) {
+	    call printf ("The graphics device is undefined\n")
+	    return
+	}
+
+	# Check for undefined input image.
+	if (im == NULL) {
+	    call printf ("The graphics device is undefined\n")
+	    return
+	}
+
+	# Check for an undefined center.
+	if (IS_INDEFR(x) || IS_INDEFR(y)) {
+	    call printf ("The surface plot center is undefined\n")
+	    return
+	}
+
+	# Get the data.
+	ny = clgeti ("surfplot.nlines")
+	nx = clgeti ("surfplot.ncolumns")
+	x1 = x - (nx - 1) / 2 + 0.5
+	x2 = x + nx / 2 + 0.5
+	y1 = y - (ny - 1) / 2 + 0.5
+	y2 = y + ny / 2 + 0.5
+	data = pt_gdata (im, x1, x2, y1, y2)
+	if (data == NULL) {
+	    call printf ("The requested image section if off the image\n")
+	    return
+	}
+
+	call smark (sp)
+
+	# Set the title.
+	call salloc (longtitle, 2 * SZ_LINE, TY_CHAR)
+	if (clgetb ("surfplot.banner")) {
+	    Memc[longtitle] = '\n'
+	    call sysid (Memc[longtitle+1], 2 * SZ_LINE)
+	    call sprintf (Memc[longtitle+strlen(Memc[longtitle])], 2 * SZ_LINE,
+	        "\nObject at x: %g  y: %g\nSurface plot of %s[%d:%d,%d:%d]")
+	        call pargr (x)
+	        call pargr (y)
+	        call pargstr (IM_HDRFILE(im))
+	        call pargi (x1)
+	        call pargi (x2)
+	        call pargi (y1)
+	        call pargi (y2)
+	} else {
+	    call sprintf (Memc[longtitle], 2 * SZ_LINE,
+	        "\nObject at x: %g  y: %g\nSurface plot of %s[%d:%d,%d:%d]")
+	        call pargr (x)
+	        call pargr (y)
+	        call pargstr (IM_HDRFILE(im))
+	        call pargi (x1)
+	        call pargi (x2)
+	        call pargi (y1)
+	        call pargi (y2)
+	}
+
+	# Initialize the plot.
+	call gclear (gd)
+
+	# Set the viewport, turn off axes drawing.
+	call gsview (gd, 0.1, 0.9, 0.1, 0.9)
+	call gseti (gd, G_DRAWAXES, NO)
+	call glabax (gd, Memc[longtitle], "", "")
+
+	nx = x2 - x1 + 1
+	ny = y2 - y1 + 1
+	npts = nx * ny
+
+	# Take floor and ceiling if enabled (nonzero).
+	floor = clgetr ("surfplot.floor")
+	ceiling = clgetr ("surfplot.ceiling")
+	if (IS_INDEFR (floor) && IS_INDEFR (ceiling))
+	    sdata = data
+	else {
+	    call salloc (sdata, npts, TY_REAL)
+	    call amovr (Memr[data], Memr[sdata], npts)
+	    if (! IS_INDEFR (floor) && ! IS_INDEFR (ceiling)) {
+		floor = min (floor, ceiling)
+		ceiling = max (floor, ceiling)
+	    }
+	    if (! IS_INDEFR (floor))
+		call amaxkr (Memr[sdata], floor, Memr[sdata], npts)
+	    if (! IS_INDEFR (ceiling))
+		call aminkr (Memr[sdata], ceiling, Memr[sdata], npts)
+	}
+
+	# Open graphics device and make plot.
+	call gopks (STDERR)
+	wkid = 1
+	call gopwk (wkid, 6, gd)
+	call gacwk (wkid)
+
+	first = 1
+	call srfabd()
+	call ggview (gd, vpx1, vpx2, vpy1, vpy2)
+	call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1)
+
+	angh = clgetr ("surfplot.angh")
+	angv = clgetr ("surfplot.angv")
+	call salloc (work, 2 * (2*nx*ny+nx+ny), TY_REAL)
+	call ezsrfc (Memr[sdata], nx, ny, angh, angv, Memr[work])
+
+	if (clgetb ("surfplot.axes")) {
+	    call gswind (gd, real (x1), real (x2), real (y1), real (y2))
+	    call gseti (gd, G_CLIP, NO)
+	    call pt_perimeter (gd, Memr[sdata], nx, ny, angh, angv)
+	}
+
+	call gdawk (wkid)
+	call gclks ()
+	call sfree (sp)
+end
+
+
+# PT_CPLOT -- Contour map
+# This is an interface to the NCAR CONREC routine.
+
+procedure pt_cplot (gd, im, x, y)
+ 
+pointer	gd		# pointer to the graphics stream
+pointer	im		# pointer to the input image
+real	x, y		# center of the contour plot
+
+int	nx, ny, x1, x2, y1, y2, nhi, dashpat, npts, ncontours, wkid, nset
+pointer	data, sp, longtitle, data1
+real	xs, xe, ys, ye, vx1, vx2, vy1, vy2
+real	zero, floor, ceiling, zmin, zmax, interval, finc
+bool	clgetb(), fp_equalr()
+int	clgeti(), btoi(), strlen()
+pointer	pt_gdata()
+real	clgetr()
+
+int	isizel, isizem, isizep, nrep, ncrt, ilab, nulbll, ioffd
+int	ioffm, isolid, nla, nlm, first
+real	xlt, ybt, side, ext, hold[5]
+common  /conflg/ first
+common  /conre4/ isizel, isizem , isizep, nrep, ncrt, ilab, nulbll, 
+         ioffd, ext, ioffm, isolid, nla, nlm, xlt, ybt, side
+common  /noaolb/ hold
+
+begin
+	# Check for undefined graphics stream.
+	if (gd == NULL) {
+	    call printf ("The graphics device is undefined\n")
+	    return
+	}
+
+	# Check for undefined input image.
+	if (im == NULL) {
+	    call printf ("The graphics device is undefined\n")
+	    return
+	}
+
+	# Check fo an undefined center.
+	if (IS_INDEFR(x) || IS_INDEFR(y)) {
+	    call printf  ("The center of the contour plot is undefined\n")
+	    return
+	}
+
+	# Get the data.
+	ny = clgeti ("cntrplot.nlines")
+	nx = clgeti ("cntrplot.ncolumns")
+	x1 = x - (nx - 1) / 2 + 0.5
+	x2 = x + nx / 2 + 0.5
+	y1 = y - (ny - 1) / 2 + 0.5
+	y2 = y + ny / 2 + 0.5
+	data = pt_gdata (im, x1, x2, y1, y2)
+	if (data == NULL) {
+	    call printf ("The image section to be contoured is off the image\n")
+	    return
+	}
+
+	call smark (sp)
+
+	# Intialize the plot
+	call gclear (gd)
+
+	# Set the WCS.
+	xs = x1
+	xe = x2
+	ys = y1
+	ye = y2
+	call gswind (gd, xs, xe, ys, ye)
+
+	nx = x2 - x1 + 1
+	ny = y2 - y1 + 1
+	if (! clgetb ("cntrplot.fill"))
+	    call gsetr (gd, G_ASPECT, real (ny-1) / real (nx-1))
+	call gseti (gd, G_ROUND, btoi (clgetb ("cntrplot.round")))
+
+	if (clgetb ("cntrplot.box")) {
+
+	    # Get number of major and minor tick marks.
+	    call gseti (gd, G_XNMAJOR, clgeti ("cntrplot.majrx"))
+	    call gseti (gd, G_XNMINOR, clgeti ("cntrplot.minrx"))
+	    call gseti (gd, G_YNMAJOR, clgeti ("cntrplot.majry"))
+	    call gseti (gd, G_YNMINOR, clgeti ("cntrplot.minry"))
+
+	    # Label tick marks on axes ?
+	    call gseti (gd, G_LABELTICKS, btoi (clgetb ("cntrplot.ticklabels")))
+
+	    # Construct the title.
+	    call salloc (longtitle, 2 * SZ_LINE, TY_CHAR)
+	    if (clgetb ("cntrplot.banner")) {
+		call sysid (Memc[longtitle], 2 * SZ_LINE)
+	        call sprintf (Memc[longtitle+strlen(Memc[longtitle])],
+		2 * SZ_LINE,
+	     "\nObject at x: %g  y: %g\n\nContour plot of %s[%d:%d,%d:%d]\n")
+		    call pargr (x)
+		    call pargr (y)
+		    call pargstr (IM_HDRFILE(im))
+		    call pargi (x1)
+		    call pargi (x2)
+		    call pargi (y1)
+		    call pargi (y2)
+	    } else {
+	        call sprintf (Memc[longtitle], 2 * SZ_LINE,
+		"\nObject at x: %g  y: %g\n\nContour plot of %s[%d:%d,%d:%d]\n")
+		    call pargr (x)
+		    call pargr (y)
+		    call pargstr (IM_HDRFILE(im))
+		    call pargi (x1)
+		    call pargi (x2)
+		    call pargi (y1)
+		    call pargi (y2)
+	    }
+
+	    call glabax (gd, Memc[longtitle], "", "")
+	}
+
+	# First of all, intialize conrec's block data before altering any
+	# parameters in common.
+	first = 1
+	call conbd
+
+	# Set the contouring parameters.
+	zero = clgetr ("cntrplot.zero")
+	floor = clgetr ("cntrplot.floor")
+	ceiling	= clgetr ("cntrplot.ceiling")
+	nhi = clgeti ("cntrplot.nhi")
+	dashpat	= clgeti ("cntrplot.dashpat")
+
+	# Resolve INDEF limits.
+	npts = nx * ny
+	if (IS_INDEFR (floor) || IS_INDEFR (ceiling)) {
+	    call alimr (Memr[data], npts, zmin, zmax)
+	    if (IS_INDEFR (floor))
+	        floor = zmin
+	    if (IS_INDEFR (ceiling))
+	        ceiling = zmax
+	}
+
+	# Apply the zero point shift.
+	if (abs (zero) > EPSILON) {
+	    call salloc (data1, npts, TY_REAL)
+	    call asubkr (Memr[data], zero, Memr[data1], npts)
+	    floor = floor - zero
+	    ceiling = ceiling - zero
+	} else
+	    data1 = data
+
+	# Avoid conrec's automatic scaling.
+	if (fp_equalr (floor, 0.0))
+	    floor = EPSILON
+	if (fp_equalr (ceiling, 0.0))
+	    ceiling = EPSILON
+
+	# The user can suppress the contour labelling by setting the common
+	# parameter "ilab" to zero.
+	if (btoi (clgetb ("cntrplot.label")) == NO)
+	    ilab = 0
+	else
+	    ilab = 1
+
+	# User can specify either the number of contours or the contour
+	# interval, or let conrec pick a nice number.  Get params and
+	# encode the FINC param expected by conrec.
+
+	ncontours = clgeti ("cntrplot.ncontours")
+	if (ncontours <= 0) {
+	    interval = clgetr ("cntrplot.interval")
+	    if (interval <= 0.0)
+		finc = 0
+	    else
+		finc = interval
+	} else
+	    finc = - abs (ncontours)
+
+	# Open device and make contour plot.
+	call gopks (STDERR)
+	wkid = 1
+	call gopwk (wkid, 6, gd)
+	call gacwk (wkid)
+
+	# Make the contour plot.
+	nset = 1	# No conrec viewport
+	ioffm = 1	# No conrec box
+	call gswind (gd, 1., real (nx), 1., real (ny))
+	call ggview (gd, vx1, vx2, vy1, vy2)
+	call set (vx1, vx2, vy1, vy2, 1.0, real (nx), 1.0, real (ny), 1)
+	call conrec (Memr[data1], nx, nx, ny, floor, ceiling, finc, nset,
+	    nhi, -dashpat)
+
+	call gdawk (wkid)
+	call gclks ()
+
+	call gswind (gd, xs, xe, ys, ye)
+	if (fp_equalr (hold[5], 1.0)) {
+	    call sprintf (Memc[longtitle], 2 * SZ_LINE, 
+                "\n\nContoured from %g to %g, interval = %g\n\n")
+	        call pargr (hold[1])
+	        call pargr (hold[2])
+	        call pargr (hold[3])
+	} else {
+	    call sprintf (Memc[longtitle], 2 * SZ_LINE,
+        "\n\nContoured from %g to %g, interval = %g, labels scaled by %g\n\n")
+	            call pargr (hold[1])
+	            call pargr (hold[2])
+	            call pargr (hold[3])
+		    call pargr (hold[5])
+	}
+
+	call gseti (gd, G_DRAWAXES, NO)
+	call glabax (gd, Memc[longtitle], "", "")
+
+	call sfree (sp)
+end
+
+
+# PT_PERIMETER -- Draw and label axes around the surface plot.
+
+procedure pt_perimeter (gd, z, ncols, nlines, angh, angv)
+
+pointer	gd			# graphics pointer
+int	ncols			# number of image columns
+int	nlines			# number of image lines
+real	z[ncols, nlines]	# array of intensity values
+real	angh			# angle of horizontal inclination
+real	angv			# angle of vertical inclination
+
+char	tlabel[10]
+int	i, j
+pointer	sp, x_val, y_val, kvec
+real	xmin, ymin, delta, fact1, flo, hi, xcen, ycen
+real	x1_perim, x2_perim, y1_perim, y2_perim, z1, z2
+real	wc1, wc2, wl1, wl2, del
+int	itoc()
+
+data  	fact1 /2.0/
+real	vpx1, vpx2, vpy1, vpy2
+common	/noaovp/ vpx1, vpx2, vpy1, vpy2
+
+begin
+	call smark (sp)
+	call salloc (x_val, ncols + 2, TY_REAL)
+	call salloc (y_val, nlines + 2, TY_REAL)
+	call salloc (kvec, max (ncols, nlines) + 2, TY_REAL)
+
+	# Get window coordinates set up calling procedure.
+	call ggwind (gd, wc1, wc2, wl1, wl2)
+
+	# Set up window, viewport for output.  The coordinates returned
+	# from trn32s are in the range [1-1024].
+	call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1)
+
+	# Find range of z for determining perspective.
+	flo = MAX_REAL
+	hi = -flo
+	do j = 1, nlines {
+	    call alimr (z[1,j], ncols, z1, z2)
+	    flo = min (flo, z1)
+	    hi = max (hi, z2)
+	}
+
+	# Set up linear endpoints and spacing as used in surface.
+        delta = (hi-flo) / (max (ncols, nlines) -1.) * fact1
+        xmin = -(real (ncols/2)  * delta + real (mod (ncols+1, 2))  * delta)
+        ymin = -(real (nlines/2) * delta + real (mod (nlines+1, 2)) * delta)
+	del = 2.0 * delta
+
+	# The perimeter is separated from the surface plot by the 
+	# width of delta.  
+	x1_perim = xmin - delta
+	y1_perim = ymin - delta
+	x2_perim = xmin + (real (ncols)  * delta)
+	y2_perim = ymin + (real (nlines) * delta)
+
+	# Set up linear arrays over full perimeter range.
+	do i = 1, ncols + 2
+	    Memr[x_val+i-1] = x1_perim + (i-1) * delta
+	do i = 1, nlines + 2
+	    Memr[y_val+i-1] = y1_perim + (i-1) * delta
+
+	# Draw and label axes and tick marks.
+	# It is important that frame has not been called after calling srface.
+	# First to draw the perimeter.  Which axes get drawn depends on the
+	# values of angh and angv.  Get angles in the range [-180, 180].
+
+	if (angh > 180.)
+	    angh = angh - 360.
+	else if (angh < -180.)
+	    angh = angh + 360.
+	if (angv > 180.)
+	    angv = angv - 360.
+	else if (angv < -180.)
+	    angv = angv + 360.
+
+	# Calculate positions for the axis labels.
+	xcen = 0.5 * (x1_perim + x2_perim)
+	ycen = 0.5 * (y1_perim + y2_perim)
+
+	if (angh >= 0.0) {
+
+	    # Case 1: xy rotation positive, looking down from above mid z.
+	    if (angv >= 0.0) {
+
+		# First draw x axis.
+		call amovkr (y2_perim, Memr[kvec], ncols + 2)
+		call pt_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1)
+		call pt_label_axis (xcen, y2_perim+del, flo, "X-AXIS", -1, -2)
+		call pt_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, 
+		    flo, ncols)
+		if (itoc (int (wc1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (xmin, y2_perim+del, flo, tlabel, -1, -2)
+		if (itoc (int (wc2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (Memr[x_val+ncols], y2_perim+del, flo, 
+		    tlabel, -1, -2)
+
+		# Now draw y axis.
+		call amovkr (x2_perim, Memr[kvec], nlines + 2)
+		call pt_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1)
+		call pt_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1)
+		call pt_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1],
+		    flo, nlines)
+		if (itoc (int (wl1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x2_perim+del, ymin, flo, tlabel, 2, -1)
+		if (itoc (int (wl2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x2_perim+del, Memr[y_val+nlines], flo, 
+		    tlabel, 2, -1)
+
+	    # Case 2: xy rotation positive, looking up from below mid z.
+	    } else {
+
+		# First draw x axis.
+		call amovkr (y1_perim, Memr[kvec], ncols + 2)
+		call pt_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1)
+		call pt_label_axis (xcen, y1_perim-del, flo, "X-AXIS", -1, 2)
+		call pt_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, 
+		    flo, ncols)
+		if (itoc (int (wc1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (xmin, y1_perim-del, flo, tlabel, -1, 2)
+		if (itoc (int (wc2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (Memr[x_val+ncols], y1_perim-del, flo, 
+		    tlabel, -1, 2)
+
+		# Now draw y axis.
+		call amovkr (x1_perim, Memr[kvec], nlines + 2)
+		call pt_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1)
+		call pt_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1)
+		call pt_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1],
+		    flo, nlines)
+		if (itoc (int (wl1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x1_perim-del, ymin, flo, tlabel, 2, 1)
+		if (itoc (int (wl2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x1_perim-del, Memr[y_val+nlines], flo, 
+		    tlabel, 2, 1)
+	    }
+	}
+
+	if (angh < 0.0) {
+
+	    # Case 3: xy rotation negative, looking down from above  mid z 
+	    # (default).
+	    if (angv > 0.0) {
+
+		# First draw x axis.
+		call amovkr (y1_perim, Memr[kvec], ncols + 2)
+		call pt_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1)
+		call pt_label_axis (xcen, y1_perim-del, flo, "X-AXIS", 1, 2)
+		call pt_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, 
+		    flo, ncols)
+		if (itoc (int (wc1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (xmin, y1_perim-del, flo, tlabel, 1, 2)
+		if (itoc (int (wc2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (Memr[x_val+ncols], y1_perim-del, flo, 
+		    tlabel, 1, 2)
+
+		# Now draw y axis.
+		call amovkr (x2_perim, Memr[kvec], nlines + 2)
+		call pt_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1)
+		call pt_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1)
+		call pt_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1],
+		    flo, nlines)
+		if (itoc (int (wl1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x2_perim+del, ymin, flo, tlabel, 2, -1)
+		if (itoc (int (wl2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x2_perim+del, Memr[y_val+nlines], flo, 
+		    tlabel, 2, -1)
+
+	    # Case 4: xy rotation negative, looking up from below mid Z.
+	    } else {
+
+		# First draw x axis.
+		call amovkr (y2_perim, Memr[kvec], ncols + 2)
+		call pt_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1)
+		call pt_label_axis (xcen, y2_perim+del, flo, "X-AXIS", 1, -2)
+		call pt_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta,
+		    flo, ncols)
+		if (itoc (int (wc1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (xmin, y2_perim+del, flo, tlabel, 1, -2)
+		if (itoc (int (wc2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (Memr[x_val+ncols], y2_perim+del, flo, 
+		    tlabel, 1, -2)
+
+		# Now draw y axis.
+		call amovkr (x1_perim, Memr[kvec], nlines + 2)
+		call pt_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1)
+		call pt_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1)
+		call pt_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1],
+		    flo, nlines)
+		if (itoc (int (wl1), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x1_perim-del, ymin, flo, tlabel, 2, 1)
+		if (itoc (int (wl2), tlabel, 10) <= 0)
+		    tlabel[1] = EOS
+		call pt_label_axis (x1_perim-del, Memr[y_val+nlines], flo, 
+		    tlabel, 2, 1)
+	    }
+	}
+
+	# Flush plotit buffer before returning.
+	call plotit (0, 0, 2)
+	call sfree (sp)
+end
+
+
+# PT_DRAW_AXIS -- Draw the axes around the plot.
+
+procedure pt_draw_axis (xvals, yvals, zval, nvals)
+
+real	xvals[nvals]
+real	yvals[nvals]
+real	zval
+int	nvals
+
+int	i
+pointer	sp, xt, yt
+real	dum
+
+begin
+	call smark (sp)
+	call salloc (xt, nvals, TY_REAL)
+	call salloc (yt, nvals, TY_REAL)
+
+	do i = 1, nvals 
+	    call trn32s (xvals[i], yvals[i], zval, Memr[xt+i-1], Memr[yt+i-1], 
+		dum, 1)
+
+	call gpl (nvals, Memr[xt], Memr[yt])
+	call sfree (sp)
+end
+
+
+# PT_LABEL_AXIS -- Label the axes.
+
+procedure pt_label_axis (xval, yval, zval, sppstr, path, up)
+
+real	xval
+real	yval
+real	zval
+char	sppstr[SZ_LINE]
+int	path
+int	up
+
+int	nchars
+int	strlen()
+%	character*64 fstr
+
+begin
+	nchars = strlen (sppstr)
+%	call f77pak (sppstr, fstr, 64)
+    	call pwrzs (xval, yval, zval, fstr, nchars, CSIZE, path, up, 0)
+end
+
+
+# PT_DRAW_TICKSX -- Draw the x tick marks.
+
+procedure pt_draw_ticksx (x, y1, y2, zval, nvals)
+
+real	x[nvals]
+real	y1, y2
+real	zval
+int	nvals
+
+int	i
+real	tkx[2], tky[2], dum
+
+begin
+	do i = 1, nvals {
+	    call trn32s (x[i], y1, zval, tkx[1], tky[1], dum, 1)
+	    call trn32s (x[i], y2, zval, tkx[2], tky[2], dum, 1)
+	    call gpl (2, tkx[1], tky[1])
+	}
+end
+
+
+# PT_DRAW_TICKSY -- Draw the y tick marks.
+
+procedure pt_draw_ticksy (x1, x2, y, zval, nvals)
+
+real	x1, x2
+real	y[nvals]
+real	zval
+int	nvals
+
+int	i
+real	tkx[2], tky[2], dum
+
+begin
+	do i = 1, nvals {
+	    call trn32s (x1, y[i], zval, tkx[1], tky[1], dum, 1)
+	    call trn32s (x2, y[i], zval, tkx[2], tky[2], dum, 1)
+	    call gpl (2, tkx[1], tky[1])
+	}
+end
+
+ 
+# PT_GDATA -- Get image data with boundary checking.
+ 
+pointer procedure pt_gdata (im, x1, x2, y1, y2)
+ 
+pointer	im			# pointer to the input image
+int	x1, x2, y1, y2		# subraster limits both input and output
+ 
+int	i, nc, nl
+pointer	imgs2r()
+errchk	imgs2r
+ 
+begin
+	nc = IM_LEN(im,1)
+	nl = IM_LEN(im,2)
+	if (IS_INDEFI (x1))
+	    x1 = 1
+	if (IS_INDEFI (x2))
+	    x2 = nc
+	if (IS_INDEFI (y1))
+	    y1 = 1
+	if (IS_INDEFI (y2))
+	    y2 = nl
+ 
+	i = max (x1, x2)
+	x1 = min (x1, x2)
+	x2 = i
+	i = max (y1, y2)
+	y1 = min (y1, y2)
+	y2 = i
+
+	if (x2 < 1 || x1 > nc || y2 < 1 || y1 > nl)
+	    return (NULL)
+	
+	x1 = max (1, x1)
+	x2 = min (nc, x2)
+	y1 = max (1, y1)
+	y2 = min (nl, y2)
+
+	return (imgs2r (im, x1, x2, y1, y2))
+end
+
+
+# PT_RADPIX -- Procedure to fetch the image pixels in an annulus around
+# a given center.
+
+int procedure pt_radpix (im, wx, wy, rin, rout, rcoords, pix)
+
+pointer	im		# pointer to IRAF image
+real	wx, wy		# center of sky annulus
+real	rin, rout	# inner and outer radius of sky annulus
+real	rcoords[ARB]	# radial coordinate array
+real	pix[ARB]	# pixel array
+
+int	i, j, ncols, nlines, c1, c2, l1, l2, npix
+pointer	buf
+real	xc1, xc2, xl1, xl2, rin2, rout2, rj2, r2
+pointer	imgs2r()
+
+begin
+	if (rout <= rin)
+	    return (0)
+
+	# Test for out of bounds sky regions.
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xc1 = wx - rout
+	xc2 = wx + rout
+	xl1 = wy - rout
+	xl2 = wy + rout
+	if (xc2 < 1.0 || xc1 > real (ncols) || xl2 < 1.0 || xl1 > real (nlines))
+	    return (0)
+
+	# Compute the column and line limits.
+	c1 = max (1.0, min (real (ncols), wx - rout)) + 0.5
+	c2 = min (real (ncols), max (1.0, wx + rout)) + 0.5
+	l1 = max (1.0, min (real (nlines), wy - rout)) + 0.5
+	l2 = min (real (nlines), max (1.0, wy + rout)) + 0.5
+
+	# Fetch the sky pixels.
+	rin2 = rin ** 2
+	rout2 = rout ** 2
+	npix = 0
+
+	do j = l1, l2 {
+	    buf = imgs2r (im, c1, c2, j, j)
+	    rj2 = (wy - j) ** 2
+	    do i = c1, c2 {
+	        r2 = (wx - i) ** 2 + rj2
+		if (r2 > rin2 && r2 <= rout2) {
+		    rcoords[npix+1] = sqrt (r2)
+		    pix[npix+1] = Memr[buf+i-c1]
+		    npix = npix + 1
+		}
+	    }
+	}
+
+	return (npix)
+end