# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<xwhen.h>
include	<config.h>
include	<imhdr.h>
include	<mach.h>
include	<error.h>
include	<ctype.h>
include	<fset.h>
include	<gset.h>
include	<mwset.h>

define	SZ_BUF		2048		# Initial pixel buffer size
define	MAX_CURVES	20		# maximum curves if overplotting
define	LIST_OP		1
define	IMAGE_OP	2


# GRAPH -- Graphing utility where input may be one or more lists (y or x,y) 
# or image sections.  Multidimensional image sections are reduced to a vector 
# by computing the projection about the indicated axis.  Many options are 
# available to personalize the plot; see the manual page for a full description.

procedure t_graph()

char	input[SZ_LINE]
pointer	x[MAX_CURVES], y[MAX_CURVES], size[MAX_CURVES]
int	npix[MAX_CURVES], ncurves

bool	append, overplot
char	device[SZ_FNAME]
int	tgrjmp[LEN_JUMPBUF], epa, old_onint, status, i

bool	clgetb()
int	fstati()
extern	tgr_onint()
common	/tgrcom/ tgrjmp

begin
	# Initialize curve pointers to NULL, in case ggplot aborts without
	# allocating any buffers.
	do i = 1, MAX_CURVES {
	    x[i] = NULL
	    y[i] = NULL
	    size[i] = NULL
	    npix[i] = NULL
	}

	if (fstati (STDIN, F_REDIR) == YES)
	    call strcpy ("STDIN", input, SZ_FNAME)
	else
	    call clgstr ("input", input, SZ_LINE)

	# Fetch plotting parameters.

	call clgstr ("device", device, SZ_FNAME)
	overplot = clgetb ("overplot")
	append = clgetb ("append")

	# Install interrupt exception handler.
	call zlocpr (tgr_onint, epa)
	call xwhen (X_INT, epa, old_onint)

	call zsvjmp (tgrjmp, status)
	if (status == OK) {
	    # Fetch remaining params and draw the plot.
	    iferr (call ggplot (device, overplot, append, input, x, y,
		size, npix, ncurves))
		status = ERR
	}

	if (status == ERR)
	    call fseti (STDOUT, F_CANCEL, OK)

	# Return buffer space whether or not an error occurs while plotting.

	do i = 1, MAX_CURVES {
	    call mfree (x[i], TY_REAL)
	    call mfree (y[i], TY_REAL)
	    call mfree (size[i], TY_REAL)
	}

	if (status == ERR)
	    call erract (EA_ERROR)
end


# TGR_ONINT -- Interrupt handler for the task graph.  Branches back to ZSVJMP
# in the main routine to permit shutdown without an error message.

procedure tgr_onint (vex, next_handler)

int	vex			# Virtual exception
int	next_handler		# not used

int	tgrjmp[LEN_JUMPBUF]
common	/tgrcom/ tgrjmp

begin
	call xer_reset()
	call zdojmp (tgrjmp, vex)
end


# GGPLOT -- Does the real work of making the graph, after the graphics
# devics has been opened.  Fetch remaining parameters, read in the data,
# and make the plot.

procedure ggplot (device, overplot, append, input, x, y, size, npix, ncurves)

char	device[SZ_FNAME]	# Graphics device
bool	overplot		# Overplot graph
bool	append			# Append graph
char	input[ARB]		# List of operands to be plotted
pointer	x[MAX_CURVES]		# X values
pointer y[MAX_CURVES]		# Y values
pointer size[MAX_CURVES]	# Size of markers to plot
int	npix[MAX_CURVES]	# Number of points per curve
int	ncurves			# Number of curves to overplot

pointer	gd
char	xlabel[SZ_LINE], ylabel[SZ_LINE], title[SZ_LINE]
char	marker[SZ_FNAME], wcs[SZ_FNAME], xformat[SZ_FNAME], yformat[SZ_FNAME]
bool	pointmode, lintran, xautoscale, yautoscale
bool	drawbox, transpose, rdmarks
int	ltype, color, ip1, ip2
int	xtran, ytran, axis, ticklabels, i, marker_type, j
real	p1, p2, q1, q2, wx1, wx2, wy1, wy2, szmarker, vx1, vx2, vy1, vy2
real	xx, yy, sz, szx, szy
pointer	sp, ltypes, colors, ptemp

pointer	gopen()
bool	clgetb(), streq(), fp_equalr()
int	clgeti(), gg_rdcurves(), ctoi(), gstati()
real	clgetr(), plt_iformatr()
errchk	clgetb, clgeti, clgstr, clgetr, glabax, gpmark
errchk	gswind, gseti, gg_rdcurves, gascale, grscale

begin
	call smark (sp)
	call salloc (ltypes, SZ_LINE, TY_CHAR)
	call salloc (colors, SZ_LINE, TY_CHAR)

	# If computing projection along an axis (collapsing a multidimensional
	# section to a vector), fetch axis number.  Get wcs string.
	axis = clgeti ("axis")
	call clgstr ("wcs", wcs, SZ_FNAME)

	# Set the line type and color lists.
	i = 0
	call clgstr ("ltypes", Memc[ltypes], SZ_LINE)
	for (ip1=ltypes; Memc[ip1]!=EOS; ip1=ip1+1) {
	    if (Memc[ip1] == ',')
		Memc[ip1] = ' '
	    if (IS_DIGIT(Memc[ip1]))
		i = i + 1
	}
	if (i == 0)
	    Memc[ltypes] = EOS
	ip1 = 1
	ltype = 0

	i = 0
	call clgstr ("colors", Memc[colors], SZ_LINE)
	for (ip2=colors; Memc[ip2]!=EOS; ip2=ip2+1) {
	    if (Memc[ip2] == ',')
		Memc[ip2] = ' '
	    if (IS_DIGIT(Memc[ip2]))
		i = i + 1
	}
	if (i == 0)
	    Memc[colors] = EOS
	ip2 = 1
	color = 0

	# If pointmode is enabled, get marker character to be used to mark
	# points.  The size of the character is given
	# by szmarker; if zero, the size will be taken from the input list.

	pointmode = clgetb ("pointmode")
	szmarker = 0.0
	rdmarks = false

	if (pointmode) {
	    call clgstr ("marker", marker, SZ_FNAME)
	    call init_marker (marker, marker_type)
	    if (marker_type != GM_POINT) {
		szmarker = clgetr ("szmarker")
		rdmarks = (szmarker <= 0)
	    }
	} else
	    call clgstr ("marker", marker, SZ_FNAME)

	# Read all the curves specified by the operands in input into memory.
	# Get the first image title and coordinate label.

	title[1] = EOS
	xlabel[1] = EOS
	ylabel[1] = EOS
	xformat[1] = EOS
	yformat[1] = EOS
	ncurves = gg_rdcurves (input, title, xlabel, ylabel, xformat,
	    x, y, size, npix, axis, wcs, rdmarks)

	if (overplot || append)
	    gd = gopen (device, APPEND, STDGRAPH)
	else
	    gd = gopen (device, NEW_FILE, STDGRAPH)

	xautoscale = false
	yautoscale = false

	# Set window and viewport.  If user window has not been set, enable
	# autoscaling.  If device viewport has not been set, let glabax
	# handle the viewport internally.

	if (!append) {
	    wx1 = clgetr ("wx1")
	    wx2 = clgetr ("wx2")
	    wy1 = clgetr ("wy1")
	    wy2 = clgetr ("wy2")

	    if (fp_equalr (wx1, wx2))
		xautoscale = true
	    if (fp_equalr (wy1, wy2))
		yautoscale = true

	    vx1 = clgetr ("vx1")
	    vx2 = clgetr ("vx2")
	    vy1 = clgetr ("vy1")
	    vy2 = clgetr ("vy2")

	    if (!(fp_equalr (vx1, vx2)) && !(fp_equalr (vy1, vy2)))
	        call gsview (gd, vx1, vx2, vy1, vy2)

	    if (!clgetb ("fill"))
	        call gseti (gd, G_ASPECT, 1)

	    if (clgetb ("round"))
	        call gseti (gd, G_ROUND, YES)
	}

	# Draw box around plot?
	drawbox = false
	if (!append)
	    if (clgetb ("box"))
		drawbox = true

	if (drawbox) {
	    # Get number of major and minor tick marks.
	    call gseti (gd, G_XNMAJOR, clgeti ("majrx"))
	    call gseti (gd, G_XNMINOR, clgeti ("minrx"))
	    call gseti (gd, G_YNMAJOR, clgeti ("majry"))
	    call gseti (gd, G_YNMINOR, clgeti ("minry"))

	    # Fetch plot title, labels and format
	    call clgstr ("title", wcs, SZ_LINE)
	    if (!streq (wcs, "imtitle"))
		call strcpy (wcs, title, SZ_LINE)

	    call clgstr ("xlabel", wcs, SZ_LINE)
	    if (!streq (wcs, "wcslabel"))
		call strcpy (wcs, xlabel, SZ_LINE)

	    call clgstr ("ylabel", ylabel, SZ_LINE)

	    call clgstr ("xformat", wcs, SZ_LINE)
	    if (!streq (wcs, "wcsformat"))
		call strcpy (wcs, xformat, SZ_FNAME)

	    call clgstr ("yformat", yformat, SZ_LINE)

	    # Label tick marks on axes?
	    ticklabels = NO
	    if (clgetb ("ticklabels"))
		ticklabels = YES
	}

	# Perform linear transformation on the X axis?
	lintran = clgetb ("lintran")
	if (lintran) {
	    p1 = clgetr ("p1")
	    p2 = clgetr ("p2")
	    q1 = clgetr ("q1")
	    q2 = clgetr ("q2")
	}

	# Transpose X,Y axes?
	transpose = clgetb ("transpose")

	# Log scale?  Call gswind to set log scaling regardless of whether
	# the user window is known; if the user window was not input,
	# autoscaling will reset it later.

	if (append) {
	    xtran = gstati (gd, G_XTRAN)
	    ytran = gstati (gd, G_YTRAN)
	    call ggwind (gd, wx1, wx2, wy1, wy2)
	} else {
	    xtran = GW_LINEAR
	    if (clgetb ("logx"))
		xtran = GW_LOG
	    ytran = GW_LINEAR
	    if (clgetb ("logy"))
		ytran = GW_LOG
	    wx1 = plt_iformatr (wx1, xformat)
	    wx2 = plt_iformatr (wx2, xformat)
	    wy1 = plt_iformatr (wy1, yformat)
	    wy2 = plt_iformatr (wy2, yformat)
	    call gswind (gd, wx1, wx2, wy1, wy2)
	    call gseti (gd, G_XTRAN, xtran)
	    call gseti (gd, G_YTRAN, ytran)
	}

	# Carry out linear transformation on X coords, if desired.
	if (lintran)
	    do i = 1, ncurves
		call gg_lintran (Memr[x[i]], npix[i], p1,p2, q1,q2)

	# Swap axes, if enabled.  Note that the linear transformation of
	# the x-axis should be performed before axes are swapped.  This is
	# because the purpose of the lintran option is to provide a means
	# of assigning a coordinate system to a pixel array.

	if (transpose)
	    do i = 1, ncurves {
		ptemp = x[i]
		x[i] = y[i]
		y[i] = ptemp
	    }

	# Autoscale if enabled.
	if (xautoscale) {
	    call gascale (gd, Memr[x[1]], npix[1], 1)
	    if (ncurves > 1) {
		do i = 2, ncurves
		    call grscale (gd, Memr[x[i]], npix[i], 1)
	    }
	} 

	if (yautoscale) {
	    call gascale (gd, Memr[y[1]], npix[1], 2)
	    if (ncurves > 1) {
		do i = 2, ncurves
		    call grscale (gd, Memr[y[i]], npix[i], 2)
	    }
	}

	# Draw box around plot if enabled.
	if (drawbox) {
	    call gsets (gd, G_XTICKFORMAT, xformat)
	    call gsets (gd, G_YTICKFORMAT, yformat)
	    call gseti (gd, G_LABELTICKS, ticklabels)
	    call glabax (gd, title, xlabel, ylabel)
	}

	# Draw the curves.
	do i = 1, ncurves {
	    if (Memc[ltypes] == EOS)
		ltype = ltype + 1
	    else if (ctoi (Memc[ltypes], ip1, j) > 0)
		ltype = j
	    ltype = mod (ltype - 1, 4) + 1
	    call gseti (gd, G_PLTYPE, ltype)
	    if (Memc[colors] == EOS)
		color = color + 1
	    else if (ctoi (Memc[colors], ip2, j) > 0)
		color = j
	    color = mod (color - 1, 9) + 1
	    call gseti (gd, G_PLCOLOR, color)
	    if (pointmode) {
		if (!rdmarks) {
		    call amovkr (szmarker, Memr[size[i]], npix[i])
		    call gpmark (gd, Memr[x[i]], Memr[y[i]], npix[i], 
		        marker_type, Memr[size[i]], Memr[size[i]])
		} else {
		    if (szmarker < 0)
		        call amulkr (Memr[size[i]], szmarker, Memr[size[i]],
			    npix[i])
		    do j = 1, npix[i] {
			xx = Memr[x[i]+j-1]
			yy = Memr[y[i]+j-1]
			sz = Memr[size[i]+j-1]
			szx= sz; szy = sz
			if (marker_type == GM_VEBAR)
			    szx = 1.0
			else if (marker_type == GM_HEBAR)
			    szy = 1.0
			call gmark (gd, xx, yy, marker_type, szx, szy)
		    }
		}
	    } else
		call hgpline (gd, Memr[x[i]], Memr[y[i]], npix[i], marker)
	}

	call gclose (gd)
	call sfree (sp)
end


# GG_RDCURVES -- Given the operand list as input, read in all the referenced
# lists and/or image sections, producing a list of vectors as output.  Return
# as the function value the number of curves.

int procedure gg_rdcurves (oplist, title, xlabel, ylabel, xformat,
	x, y, size, npix, axis, wcs, rdmarks)

char	oplist[ARB]		# Operand list
char	title[ARB]		# Title
char	xlabel[ARB]		# X label
char	ylabel[ARB]		# Y label
char	xformat[ARB]		# WCS coordinate format
pointer	x[ARB]			# Pointer to x vector
pointer	y[ARB]			# Pointer to y vector
pointer	size[ARB]		# Pointer to vector of marker sizes
int	npix[ARB]		# Number of values per vector
int	axis			# Axis for projection
char	wcs[ARB]		# WCS type
bool	rdmarks			# Read marks from list?

char	operand[SZ_FNAME]
int	ncurves, i, fd
int	gg_rdcurve(), imtopen(), imtgetim()

begin
	ncurves = 0

	# Read all the curves into memory.

	fd = imtopen (oplist)
	while (imtgetim (fd, operand, SZ_FNAME) != EOF) {
	    ncurves = ncurves + 1
	    if (ncurves > MAX_CURVES)
		call error (0, "Maximum of 20 curves can be overplotted")
	    i = ncurves
	    iferr {
		npix[i] = gg_rdcurve (operand, title, xlabel, ylabel,
		    xformat, x[i], y[i], size[i], axis, wcs, rdmarks)
	    } then {
		call erract (EA_WARN)
		ncurves = ncurves - 1
	    }
	}

	call imtclose (fd)

	if (ncurves == 0)
	    call error (0, "No curves read")
	else
	    return (ncurves)
end


# GG_RDCURVE -- Read a curve into memory.  The operand may specify either
# list or image input; we determine which and then call the appropriate
# input routine to access the data.  Set the image title and coordinate
# label if not previously defined.

int procedure gg_rdcurve (operand, title, xlabel, ylabel, xformat,
	x, y, size, axis, wcs, rdmarks)

char	operand[ARB]		# List of operaands to be plotted
char	title[ARB]		# Title
char	xlabel[ARB]		# X label
char	ylabel[ARB]		# Y label
char	xformat[ARB]		# WCS coordinate format
pointer	x, y, size		# Pointers to x, y and size arrays
int	axis			# Axis of image projection
char	wcs[ARB]		# WCS type
bool	rdmarks			# Read marks from list?

int	gg_rdlist2(), gg_rdimage2(), gg_optype()
errchk	gg_rdlist2, gg_rdimage2, gg_optype

begin
	if (gg_optype (operand) == LIST_OP)
	    return (gg_rdlist2 (operand, x, y, size, rdmarks))
	else
	    return (gg_rdimage2 (operand, title, xlabel, ylabel, xformat,
		x, y, size, axis, wcs))
end


# GG_OPTYPE -- Determine whether the operand argument is an image section
# or a list.  If the string is STDIN, it is a list; if a subscript is
# present, it is an image; otherwise we must test whether or not it is a
# binary file and make the decision based on that.

int procedure gg_optype (operand)

char	operand[ARB]		# Operand to be plotted
int	first, last, ip
int	access(), strncmp()

begin
	# Strip off any whitespace at the beginning or end of the string.
	for (ip=1;  IS_WHITE(operand[ip]);  ip=ip+1)
	    ;
	first = ip
	for (last=ip;  operand[ip] != EOS;  ip=ip+1)
	    if (!IS_WHITE(operand[ip]))
		last = ip

	if (first == last)
	    return (LIST_OP)
	else if (strncmp (operand[first], "STDIN", 5) == 0)
	    return (LIST_OP)
	else if (operand[last] == ']')
	    return (IMAGE_OP)
	else if (access (operand, 0, TEXT_FILE) == YES)
	    return (LIST_OP)
	else 
	    return (IMAGE_OP)
end


# GG_RDIMAGE2 -- Read an image section and compute the projection about
# one dimension, producing x and y vectors as output.  Set the title
# and coordinate label if not previously defined.

int procedure gg_rdimage2 (imsect, title, xlabel, ylabel, xformat, x, y, size,
	axis, wcs)

char	imsect[ARB]		# Image section to be plotted
char	title[ARB]		# Image title
char	xlabel[ARB]		# Coordinate label
char	ylabel[ARB]		# Pixel value label
char	xformat[ARB]		# WCS coordinate format
pointer	x, y, size		# Pointer to x, y and size vector
int	axis			# Axis about which the projection is to be taken
char	wcs[ARB]		# WCS type

int	npix, i, stridxs()
pointer	sp, im, mw, ct, axvals, str
pointer	immap(), mw_openim(), mw_sctran()
errchk	immap, im_projection, malloc, mw_openim, mw_sctran, plt_wcs

begin
	call smark (sp)
	call salloc (axvals, IM_MAXDIM, TY_REAL)
	call salloc (str, SZ_FNAME, TY_CHAR)

	im = immap (imsect, READ_ONLY, 0)

	if (axis < 1 || axis > IM_NDIM(im))
	    call error (2, "Attempt to take projection over nonexistent axis")
	npix = IM_LEN(im,axis)

	call malloc (y, npix, TY_REAL)
	call im_projection (im, Memr[y], npix, axis)

	iferr {
	    call malloc (x, npix, TY_REAL)
	    call malloc (size, npix, TY_REAL)
	} then
	    call erract (EA_FATAL)

	# Set title if not previously defined
	if (title[1] == EOS) {
	    call strcpy (IM_TITLE(im), title, SZ_LINE)
	    if (stridxs ("\n", title) == 0)
		call strcat ("\n", title, SZ_LINE)
	    call imgsection (imsect, Memc[str], SZ_LINE)
	    if (Memc[str] != EOS)
		    call strcat (Memc[str], title, SZ_LINE)
	}

	# Set WCS coordinates
	mw = mw_openim (im)
	call mw_seti (mw, MW_USEAXMAP, NO)
	ct = mw_sctran (mw, "logical", wcs, 0)
	call strcpy (wcs, Memc[str], SZ_LINE)
	do i = 1, IM_NDIM(im)
	    Memr[axvals+i-1] = (1 + IM_LEN(im, i)) / 2.
	call plt_wcs (im, mw, ct, axis, Memr[axvals], 1., real(npix), Memr[x],
	    npix, Memc[str], xformat,  SZ_FNAME)
	if (xlabel[1] == EOS)
	    call strcpy (Memc[str], xlabel, SZ_LINE)
	call mw_close (mw)

	call imunmap (im)

	call sfree (sp)
	return (npix)
end


# GG_RDLIST2 -- Read a list of two dimensional data pairs into two type
# real arrays in memory.  Return pointers to the arrays and a count of the
# number of pixels.  If mark sizes are to be read from the input list,
# a third array of mark sizes is returned.  Mark sizes can only be given
# in two column (x,y) mode, with the mark size given as a third column.

int procedure gg_rdlist2 (fname, x, y, size, rdmarks)

char	fname[ARB]		# Name of list file
pointer	x, y, size		# Pointers to x, y and size vectors
bool	rdmarks			# Read markers from file?

int	buflen, n, fd, ncols, lineno
pointer	sp, lbuf, ip
real	xval, yval, szmark
int	getline(), nscan(), open()
errchk	open, sscan, getline, malloc

begin
	call smark (sp)
	call salloc (lbuf, SZ_LINE, TY_CHAR)

	fd = open (fname, READ_ONLY, TEXT_FILE)

	buflen = SZ_BUF
	iferr {
	    call malloc (x, buflen, TY_REAL)
	    call malloc (y, buflen, TY_REAL)
	    call malloc (size, buflen, TY_REAL)
	} then
	    call erract (EA_FATAL)

	n = 0
	ncols = 0
	lineno = 0
	szmark = 1E-2

	while (getline (fd, Memc[lbuf]) != EOF) {
	    # Skip comment lines and blank lines.
	    lineno = lineno + 1
	    if (Memc[lbuf] == '#')
		next
	    for (ip=lbuf;  IS_WHITE(Memc[ip]);  ip=ip+1)
		;
	    if (Memc[ip] == '\n' || Memc[ip] == EOS)
		next

	    # Decode the points to be plotted.
	    call sscan (Memc[ip])
		call gargr (xval)
		call gargr (yval)
		if (rdmarks)
		    call gargr (szmark)

	    # The first line determines whether we have an x,y list or a
	    # y-list.  It is an error if only one value can be decoded when
	    # processing a two column list.

	    if (ncols == 0 && nscan() > 0)
		ncols = nscan()
	    
	    switch (nscan()) {
	    case 0:
		call eprintf ("no args; %s, line %d: %s\n")
		    call pargstr (fname)
		    call pargi (lineno)
		    call pargstr (Memc[lbuf])
		next
	    case 1:
		if (ncols == 2) {
		    call eprintf ("only 1 arg; %s, line %d: %s\n")
			call pargstr (fname)
			call pargi (lineno)
			call pargstr (Memc[lbuf])
		    next
		} else {
		    yval = xval
		    xval = n + 1.0
		}
	    case 2:
		if (rdmarks) {
		    call eprintf ("no szmark field; %s, line %d: %s\n")
			call pargstr (fname)
			call pargi (lineno)
			call pargstr (Memc[lbuf])
		    szmark = 1E-2
		}
	    }

	    n = n + 1
	    if (n > buflen) {
		buflen = buflen + SZ_BUF
		call realloc (x, buflen, TY_REAL)
		call realloc (y, buflen, TY_REAL)
		call realloc (size, buflen, TY_REAL)
	    }

	    Memr[x+n-1] = xval
	    Memr[y+n-1] = yval
	    if (rdmarks)
		Memr[size+n-1] = szmark
	}

	call realloc (x, n, TY_REAL)
	call realloc (y, n, TY_REAL)
	call realloc (size, n, TY_REAL)

	call close (fd)
	call sfree (sp)
	return (n)
end


# GG_LINTRAN -- Linearly transform a vector.  Map pixel values P1,P2
# onto Q1,Q2.

procedure gg_lintran (x, npix, p1in, p2in, q1, q2)

real	x[npix]			# Vector to transform
int	npix			# Number of pixels in vector
real	p1in, p2in		# Range of input values to map
real	q1, q2			# Range for output values
real	p1, p2
real	xscale

begin
	# If P1 and P2 are not set, use full range of input pixels indices.
	if (p1in == 0 && p2in == 0) {
	    p1 = 1.0
	    p2 = npix
	} else {
	    p1 = p1in
	    p2 = p2in
	}

	if (p2 - p1 == 0)
	    xscale = (q2 - q1)
	else
	    xscale = (q2 - q1) / (p2 - p1)

	call asubkr (x, p1, x, npix)
	call amulkr (x, xscale, x, npix)
	call aaddkr (x, q1, x, npix)
end