Repatch (from linux) of OSX IRAF

1 files changed, 501 insertions, 0 deletions

diff --git a/pkg/plot/t_surface.x b/pkg/plot/t_surface.x
new file mode 100644
index 00000000..7faa38a2
--- /dev/null
+++ b/pkg/plot/t_surface.x
@@ -0,0 +1,501 @@
+include	<error.h>
+include	<mach.h>
+include	<gset.h>
+include	<config.h>
+include	<xwhen.h>
+include	<imhdr.h>
+include	<fset.h>
+
+define	DUMMY		6
+define	SZ_TLABEL	10
+define	CSIZE		24
+
+
+# SURFACE -- Draw a perspective view of an image section.  The altitude
+# and azimuth of the viewing angle are variable.  Floor and ceiling
+# constraints may be applied to the image data before plotting if desired.
+
+procedure t_surface()
+
+char	imsect[SZ_FNAME]
+char	device[SZ_FNAME], title[SZ_LINE]
+bool	label, sub, pre
+pointer	im, subras, work
+int	ncols, nlines, mode, wkid, nx, ny, npix
+int	epa, status, old_onint, tsujmp[LEN_JUMPBUF]
+int	xres, yres, first
+real	angh, angv, imcols, imlines
+real	floor, ceiling, vpx1, vpx2, vpy1, vpy2
+
+pointer gp, gopen()
+bool	clgetb(), streq()
+int	clgeti(), surf_limits()
+real	clgetr()
+extern	tsu_onint()
+pointer	immap(), plt_getdata()
+common	/tsucom/ tsujmp
+common  /noaovp/ vpx1, vpx2, vpy1, vpy2
+common  /frstfg/ first
+
+define	exit_ 91
+begin
+	# First initialize srface common blocks before changing any parameters
+	first = 1
+	call srfabd
+
+	# Get image section string and output device.
+	call clgstr ("image", imsect, SZ_FNAME)
+	call clgstr ("device", device, SZ_FNAME)
+
+	# Map image and open graphics device.
+	im = immap (imsect, READ_ONLY, 0)
+
+	angh = clgetr ("angh")
+	angv = clgetr ("angv")
+	floor = clgetr ("floor")
+	ceiling = clgetr ("ceiling")
+	label = clgetb ("label")
+
+	call clgstr ("title", title, SZ_LINE)
+	if (streq (title, "imtitle")) {
+	    call strcpy (imsect, title, SZ_LINE)
+	    call strcat (": ", title, SZ_LINE)
+	    call strcat (IM_TITLE(im), title, SZ_LINE)
+	}
+
+	# If a label is to be drawn, don't use the full device viewport for
+	# the surface plot.  This doesn't allow room for the axes and labels.
+
+	if (label) {
+	    vpx1 = 0.10
+	    vpx2 = 0.90
+	    vpy1 = 0.10
+	    vpy2 = 0.90
+	} 
+
+	mode = NEW_FILE
+	if (clgetb ("append"))
+	    mode = APPEND
+
+	ncols  = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	imcols = real (ncols)
+	imlines = real (nlines)
+
+	xres = clgeti ("xres")
+	yres = clgeti ("yres")
+	sub = clgetb ("subsample")
+	pre = clgetb ("preserve")
+	
+	# Get data with proper resolution.  Procedure plt_getdata returns a
+	# pointer to the data matrix to be contoured.  The resolution is
+	# decreased by the specified method in this procedure.  The image
+	# header pointer can be unmapped after plt_getdata is called.
+
+	nx = 0
+	ny = 0
+	subras = plt_getdata (im, sub, pre, xres, yres, nx, ny)
+	call imunmap (im)
+
+	# Allocate the working storage needed by EZSRFC.
+	#call malloc (work, (2 * nx * ny) + nx + ny, TY_REAL)
+	call malloc (work, 2 * ((2 * nx * ny) + nx + ny), TY_REAL)
+
+	# Take floor and ceiling if enabled (nonzero).
+	npix = nx * ny
+	if (surf_limits (Memr[subras], npix, floor, ceiling) == ERR)
+	    goto exit_
+
+	# Open graphics device and make plot.
+	call gopks (STDERR)
+	wkid = 1
+	gp = gopen (device, mode, STDGRAPH)
+	call gopwk (wkid, DUMMY, gp)
+	call gacwk (wkid)
+	call gtext (gp, 0.5, .96, title, "s=0.8;f=b;h=c")
+	call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1)
+
+	# Install interrupt exception handler.
+	call zlocpr (tsu_onint, epa)
+	call xwhen (X_INT, epa, old_onint)
+
+	call zsvjmp (tsujmp, status)
+	if (status == OK)
+	    call ezsrfc (Memr[subras], nx, ny, angh, angv, Memr[work])
+	else {
+	    call gcancel (gp)
+	    call fseti (STDOUT, F_CANCEL, OK)
+	}
+
+	if (label) {
+	    # Establish plotting window in full scale image coordinates.
+	    call gswind (gp, 1.0, imcols, 1.0, imlines)
+	    call gseti (gp, G_CLIP, NO)
+	    call srf_perimeter (gp, Memr[subras], nx, ny, angh, angv)
+	}
+
+	call gdawk (wkid)
+	call gclwk (wkid)
+	call gclks ()
+exit_
+	call mfree (subras, TY_REAL)
+	call mfree (work, TY_REAL)
+
+end
+
+
+# TSU_ONINT -- Interrupt handler for the task surface.  Branches back to ZSVJMP
+# in the main routine to permit shutdown without an error message.
+
+procedure tsu_onint (vex, next_handler)
+
+int	vex		# virtual exception
+int	next_handler	# not used
+
+int	tsujmp[LEN_JUMPBUF]
+common	/tsucom/ tsujmp
+
+begin
+	call xer_reset()
+	call zdojmp (tsujmp, vex)
+end
+
+
+# SURF_LIMITS -- Apply the floor and ceiling constraints to the subraster.
+# If both values are exactly zero, they are not applied.
+
+int procedure surf_limits (ras, npix, floor, ceiling)
+
+real	ras[npix]		# Input array of pixels
+int	npix			# npixels in array
+real	floor, ceiling		# user specified parameters
+int	apply, i
+real	tfloor, tceiling, delta, dmin, dmax
+bool	fp_equalr()
+
+begin
+	tfloor = floor
+	tceiling = ceiling
+	apply = YES
+
+	call alimr (ras, npix, dmin, dmax)
+	if (fp_equalr (dmin, dmax)) {
+	    call eprintf ("Constant valued array; no plot drawn\n")
+	    return (ERR)
+	}
+
+	if (fp_equalr (tfloor, INDEF))
+	    tfloor = dmin
+	if (fp_equalr (tceiling, INDEF))
+	    tceiling = dmax
+
+	delta = tceiling - tfloor
+	if (delta < 0.0) {
+	    # specified ceiling is lower than floor, flip them
+	    floor = tceiling
+	    ceiling = tfloor
+	} else if (fp_equalr (delta, 0.0)) {
+	    # degenerate values
+	    apply = NO
+	    floor = dmin
+	    ceiling = dmax
+	    call eprintf (
+ "Floor and ceiling are degenerate values and will be ignored\n")
+	} else {
+	    # Non-degenerate, ceiling exceedes floor as expected
+	    floor = tfloor
+	    ceiling = tceiling
+	}
+
+	if (apply == YES) {
+	    # First verify that floor and ceiling are valid
+	    if (dmax <= floor) {
+	        call eprintf ("Entire image is at or below specified floor\n")
+		return (ERR)
+	    }
+	    if (dmin >= ceiling) {
+	        call eprintf ("Entire image is at or above specified ceiling\n")
+		return (ERR)
+	    }
+
+	    do i = 1, npix {
+		# Apply surface limits
+		ras[i] = max (floor, ras[i])
+		ras[i] = min (ceiling, ras[i])
+	    }
+	}
+	 return (OK)
+end
+
+
+# SRF_PERIMETER -- draw and label axes around the surface plot.
+
+procedure srf_perimeter (gp, z, ncols, nlines, angh, angv)
+
+pointer	gp			# 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
+
+pointer	sp, x_val, y_val, kvec
+char	tlabel[SZ_TLABEL]
+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	i, j, junk
+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 in calling procedure.
+	call ggwind (gp, 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) {
+	    if (angv >= 0) {
+		# Case 1: xy rotation positive, looking down from above mid Z
+
+		# First draw x axis
+		call amovkr (y2_perim, Memr[kvec], ncols + 2)
+		call srf_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1)
+		call srf_label_axis (xcen, y2_perim+del, flo, "X-AXIS", -1, -2)
+		call srf_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, 
+		    flo, ncols)
+		call srf_label_axis (xmin, y2_perim+del, flo, "1", -1, -2)
+		junk = itoc (int (wc2), tlabel, SZ_TLABEL)
+		call srf_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 srf_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1)
+		call srf_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1)
+		call srf_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1],
+		    flo, nlines)
+		call srf_label_axis (x2_perim+del, ymin, flo, "1", 2, -1)
+		junk = itoc (int (wl2), tlabel, SZ_TLABEL)
+		call srf_label_axis (x2_perim+del, Memr[y_val+nlines], flo, 
+		    tlabel, 2, -1)
+	    } else {
+		# Case 2: xy rotation positive, looking up from below mid Z
+		# First draw x axis
+		call amovkr (y1_perim, Memr[kvec], ncols + 2)
+		call srf_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1)
+		call srf_label_axis (xcen, y1_perim-del, flo, "X-AXIS", -1, 2)
+		call srf_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, 
+		    flo, ncols)
+		call srf_label_axis (xmin, y1_perim-del, flo, "1", -1, 2)
+		junk = itoc (int (wc2), tlabel, SZ_TLABEL)
+		call srf_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 srf_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1)
+		call srf_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1)
+		call srf_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1],
+		    flo, nlines)
+		call srf_label_axis (x1_perim-del, ymin, flo, "1", 2, 1)
+		junk = itoc (int (wl2), tlabel, SZ_TLABEL)
+		call srf_label_axis (x1_perim-del, Memr[y_val+nlines], flo, 
+		    tlabel, 2, 1)
+	    }
+	}
+
+	if (angh < 0) {
+	    if (angv > 0) {
+		# Case 3: xy rotation negative, looking down from above  mid Z 
+		# (default).  First draw x axis
+		call amovkr (y1_perim, Memr[kvec], ncols + 2)
+		call srf_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1)
+		call srf_label_axis (xcen, y1_perim-del, flo, "X-AXIS", 1, 2)
+		call srf_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, 
+		    flo, ncols)
+		call srf_label_axis (xmin, y1_perim-del, flo, "1", 1, 2)
+		junk = itoc (int (wc2), tlabel, SZ_TLABEL)
+		call srf_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 srf_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1)
+		call srf_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1)
+		call srf_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1],
+		    flo, nlines)
+		call srf_label_axis (x2_perim+del, ymin, flo, "1", 2, -1)
+		junk = itoc (int (wl2), tlabel, SZ_TLABEL)
+		call srf_label_axis (x2_perim+del, Memr[y_val+nlines], flo, 
+		    tlabel, 2, -1)
+	    } else {
+		# Case 4: xy rotation negative, looking up from below mid Z
+		# First draw x axis
+		call amovkr (y2_perim, Memr[kvec], ncols + 2)
+		call srf_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1)
+		call srf_label_axis (xcen, y2_perim+del, flo, "X-AXIS", 1, -2)
+		call srf_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta,
+		    flo, ncols)
+		call srf_label_axis (xmin, y2_perim+del, flo, "1", 1, -2)
+		junk = itoc (int (wc2), tlabel, SZ_TLABEL)
+		call srf_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 srf_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1)
+		call srf_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1)
+		call srf_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1],
+		    flo, nlines)
+		call srf_label_axis (x1_perim-del, ymin, flo, "1", 2, 1)
+		junk = itoc (int (wl2), tlabel, SZ_TLABEL)
+		call srf_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
+
+
+procedure srf_draw_axis (xvals, yvals, zval, nvals)
+
+int	nvals
+real	xvals[nvals]
+real	yvals[nvals]
+real	zval
+pointer	sp, xt, yt
+int	i
+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
+
+
+procedure srf_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
+
+
+procedure srf_draw_ticksx (x, y1, y2, zval, nvals)
+
+int	nvals
+real	x[nvals]
+real	y1, y2
+real	zval
+
+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
+
+
+procedure srf_draw_ticksy (x1, x2, y, zval, nvals)
+
+int	nvals
+real	x1, x2
+real	y[nvals]
+real	zval
+
+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