"ldhplot"<-
function(data = NULL, pcoef = NULL, device = NULL)
{
# CONTROL OF PLOTTING: PRIMARY, SECONDARY, LACTATE INHIBITION
#
# data = matrix of data
#	columns:
#		"v0"     "v0calc" "weight" "[A]"    "[B]"    "[P]"
#	rows:
#		values for each data point; 
#		P = lactate, B = pyruvate, A = NADH
#		should have entries for all columns, owing to setup of 
#			ldhplot.pri.A, etc.;
#		a column of values, as for "[P]", may be all 0's;
#		
# pcoef = one-row matrix of paramter values
#	columns:
# "Vmax"     "KmA"      "KmB"      "KmAB"     "KmQ/KmPQ" "1/KiP"    "KIB_K3K4"
#	one row:
#		value for each parameter of rate law, with product inhibition
#			by P, and with abortive AP and BQ complexes;
#		should have entries for all columns, owing to setup of 
#			ldhplot.pri.A, etc.;
#		parameters 5, 6, and 7 may be made negligible by setting them
#			at 1.e-10 (do not set at zero, to avoid
#			divide-by-zero)
#	example of pcoef:
# > zz.p
#           Vmax       KmA       KmB      KmAB KmQ/KmPQ 1/KiP KIB_K3K4
# [1,] 1.222e-06 9.545e-06 0.0007459 2.523e-09    1e-10 1e-10    1e-10
#
#
# device = type of output
#	default == NULL:	X11 plot window
#	"postscript":		postscript plot to printer (currently, lp)
#	"file":			postscript plot stored as "job_pid.ps";
#				use "file" for hardcopy, so can review before
#				printing
# NOTE: if run the routines ldhplot.pri.A separately, then the
# paramters are as above for data and pcoef, but device must be, if not
# the NULL default, then the full postscript function, e.g.:
#	postscript(file="", command="lp")
#	postscript(file="", command="cat >$$.ps")
#	postscript(file="file_name")
# or sone other graphics driver (fig, iris, etc.)
#
# NOTE: defaults for data and pcoef are elements of a list, 
# 	aa$data and aa$coef
#
#
	if(is.null(data)) data <- aa$data
	if(is.null(pcoef))
		pcoef <- aa$coef
	#graphics.off()
	device.out <- device
	if(!is.null(device))
		if(device == "postscript")
			postscript(file="", command = "lp", horizontal = F)
		else if(device == "file")
			postscript(file="", command = "cat >ldhplot.pri.B-$$.ps", horizontal = F)
	ldhplot.pri.B(data, pcoef, device = device.out)
	cat("next plot? ")
	ans <- readline()
	if(ans == "n" | ans == "no" | ans == "q" | ans == "quit")
		return()
	graphics.off()
	if(!is.null(device))
		if(device == "postscript")
			postscript(file="", command = "lp", horizontal = F)
		else if(device == "file")
			postscript(file="", command = "cat >ldhplot.pri.A-$$.ps", horizontal = F)
	ldhplot.pri.A(data, pcoef, device = device.out)
	cat("next plot? ")
	ans <- readline()
	if(ans == "n" | ans == "no" | ans == "q" | ans == "quit")
		return()
	graphics.off()
	if(!is.null(device))
		if(device == "postscript")
			postscript(file="", command = "lp", horizontal = T)
		else if(device == "file")
			postscript(file="", command = "cat >ldhplot.sec-$$.ps", horizontal = T)
	ldhplot.sec(data, pcoef, device = device.out)
	cat("next plot? ")
	ans <- readline()
	if(ans == "n" | ans == "no" | ans == "q" | ans == "quit")
		return()
	graphics.off()
	if(!is.null(device))
		if(device == "postscript")
			postscript(file="", command = "lp", horizontal = T)
		else if(device == "file")
			postscript(file="", command = "cat >ldhplot.lac-$$.ps", horizontal = T)
	ldhplot.lac(data, pcoef, device = device.out)
	if(!is.null(device))
		if(device == "postscript")
			graphics.off()
		else if(device == "file")
			graphics.off()
	invisible()
}
"ldhplot.lac"<-
function(data = NULL, pcoef = NULL, device = NULL)
{
# extract data for +- lactate product inhibition
# == all data for highest concentration of A, variable B, P = {0, 150}
#
# set up plot to include axes 0, with dummy (0,0) point, but plot no points
# then plot data points
# then calculate y-axis intercepts +- P, and draw theoretical lines
#
	maxA <- max(data[, "[A]"])
	maxB <- max(data[, "[B]"])
	minB <- min(data[, "[B]"])
	maxP <- max(data[, "[P]"])
	data.P <- data[data[, "[A]"] == maxA,  ]
	xx <- 1/data.P[, "[B]"]
	yy <- 1/data.P[, "v0"]
	if(is.null(device))
		x11(width=9.2,height=6.4)	
		#x11(geometry = "920x640", close = T)	# sgi screen
	else device
	plot(c(0, xx), c(0, yy), type = "n", xlab = "", ylab = "", las = 1)
	axis(2, labels = F, tick = T, pos = 0)
	axis(1, labels = F, tick = T, pos = 0)
	title(las = 0, ylab = "1/v0", xlab = "1/[B]")
	#axes(axes = F, las = 0, ylab = "1/v0")
	#axes(axes = F, las = 0, xlab = "1/[B]")
	points(xx, yy, pch = "+")
	yynoP <- 1/pcoef[, "Vmax"] * (1 + pcoef[, "KmA"]/maxA)
	yyP <- 1/pcoef[, "Vmax"] * (1 + pcoef[, "KmA"]/maxA + pcoef[, 
		"KmQ/KmPQ"] * maxP)
	yyminB <- 1/data.P[data.P[, "[B]"] == minB, "v0calc"]
	lines(c(0, 1/minB), c(yynoP, yyminB[1]))
	lines(c(0, 1/minB), c(yyP, yyminB[length(yyminB)]))
	invisible()
}
"ldhplot.pri.A"<-
function(data = NULL, pcoef = NULL, device = NULL)
{
# PRIMARY PLOTS
# calculate points of intersection for 1/B and 1/A plots
# select only points from data array with no P == product inhibitor
# for later use, find smallest concentrations of A and B, which are the
# 	highest points of the reciprocal plots
# also for later use, extract array with replicates eliminated
# open device, after closing whatever may be open
# 1/v0 - 1/[B] plots (or same for 1/v0 - 1/[A] plots
# set up plot axes by use of theory values == v0calc, 
# 	but don't plot the points (type="n")
# add interior axes, along zeros of 1/v0, 1/concentration
# plot data points for variable B
# draw theoretical lines, from intersection to highest 1/B value for
# 	each concentration of fixed A; this involves looping over the set
# 	of fixed concentrations of A; the values are identified in the 
# 	array of uniq data points, by their sequence number
# pause between plots
#
	B.int <-  - (pcoef[, "KmA"]/pcoef[, "KmAB"])
	A.int <-  - (pcoef[, "KmB"]/pcoef[, "KmAB"])
	v.int <- 1/pcoef[, "Vmax"] * (1 - ((pcoef[, "KmA"] * pcoef[, "KmB"])/
		pcoef[, "KmAB"]))
	bb <- data[data[, "[P]"] == 0,  ]
	minB <- min(bb[, "[B]"])
	minA <- min(bb[, "[A]"])
	bb.uniq <- bb[!duplicated(paste(bb[, "[A]"], bb[, "[B]"], bb[, "[P]"])),
		]
	if(is.null(device))
		x11(width=6.4,height=9.2)	
		#x11(geometry = "640x920", close = T)	# sgi screen
	else device
	xx <- c(A.int, 1/bb.uniq[, "[A]"])
	yy <- c(v.int, 1/bb.uniq[, "v0calc"])
	plot(xx, yy, las = 1, type = "n", xlab = "", ylab = "")
	axis(2, labels = F, tick = T, pos = 0)
	axis(1, labels = F, tick = T, pos = 0)
	title(las = 0, ylab = "1/v0", xlab = "1/[A]")
	#axes(axes = F, las = 0, ylab = "1/v0")
	#axes(axes = F, las = 0, xlab = "1/[A]")
	points(1/bb[, c("[A]", "v0")], pch = "+")
	Aseq <- seq(along = bb.uniq[, 1])[bb.uniq[, "[A]"] == minA]
	for(i in Aseq)
		lines(c(A.int, 1/bb.uniq[i, "[A]"]), c(v.int, 1/bb.uniq[i, 
			"v0calc"]))
	invisible()
}
"ldhplot.pri.B"<-
function(data = NULL, pcoef = NULL, device = NULL)
{
# PRIMARY PLOTS
# calculate points of intersection for 1/B and 1/A plots
# select only points from data array with no P == product inhibitor
# for later use, find smallest concentrations of A and B, which are the
# 	highest points of the reciprocal plots
# also for later use, extract array with replicates eliminated
# open device, after closing whatever may be open
# 1/v0 - 1/[B] plots (or same for 1/v0 - 1/[A] plots
# set up plot axes by use of theory values == v0calc, 
# 	but don't plot the points (type="n")
# add interior axes, along zeros of 1/v0, 1/concentration
# plot data points for variable B
# draw theoretical lines, from intersection to highest 1/B value for
# 	each concentration of fixed A; this involves looping over the set
# 	of fixed concentrations of A; the values are identified in the 
# 	array of uniq data points, by their sequence number
#
	B.int <-  - (pcoef[, "KmA"]/pcoef[, "KmAB"])
	A.int <-  - (pcoef[, "KmB"]/pcoef[, "KmAB"])
	v.int <- 1/pcoef[, "Vmax"] * (1 - ((pcoef[, "KmA"] * pcoef[, "KmB"])/
		pcoef[, "KmAB"]))
	bb <- data[data[, "[P]"] == 0,  ]
	minB <- min(bb[, "[B]"])
	minA <- min(bb[, "[A]"])
	bb.uniq <- bb[!duplicated(paste(bb[, "[A]"], bb[, "[B]"], bb[, "[P]"])),
		]
	if(is.null(device))
		x11(width=6.4,height=9.2)	
		#x11(geometry = "640x920", close = T)	# sgi screen
	else device
	xx <- c(B.int, 1/bb.uniq[, "[B]"])
	yy <- c(v.int, 1/bb.uniq[, "v0calc"])
	plot(xx, yy, las = 1, type = "n", xlab = "", ylab = "")
	axis(2, labels = F, tick = T, pos = 0)
	axis(1, labels = F, tick = T, pos = 0)
	title(las = 0, ylab = "1/v0", xlab = "1/[B]")
	#axes(axes = F, las = 0, ylab = "1/v0")
	#axes(axes = F, las = 0, xlab = "1/[B]")
	points(1/bb[, c("[B]", "v0")], pch = "+")
	Bseq <- seq(along = bb.uniq[, 1])[bb.uniq[, "[B]"] == minB]
	for(i in Bseq)
		lines(c(B.int, 1/bb.uniq[i, "[B]"]), c(v.int, 1/bb.uniq[i, 
			"v0calc"]))
	invisible()
}
"ldhplot.sec"<-
function(data = NULL, pcoef = NULL, device = NULL)
{
# SECONDARY PLOTS
#
	B.int <- -1/pcoef[, "KmB"]
	A.int <- -1/pcoef[, "KmA"]
	bb <- data[data[, "[P]"] == 0,  ]
	minB <- min(bb[, "[B]"])
	minA <- min(bb[, "[A]"])
	B.uniq <- bb[!duplicated(paste(bb[, "[B]"])),  ]
	A.uniq <- bb[!duplicated(paste(bb[, "[A]"])),  ]
	bb.uniq <- bb[!duplicated(paste(bb[, "[A]"], bb[, "[B]"], bb[, "[P]"])),
		]
	if(is.null(device))
		x11(width=9.2,height=6.4)	
		#x11(geometry = "920x640", close = T)	#sgi screen
	else device
	xx.B <- 1/B.uniq[, "[B]"]
	yy.B <- 1/pcoef[, "Vmax"] * (1 + pcoef[, "KmB"]/B.uniq[, "[B]"])
	xx.A <- 1/A.uniq[, "[A]"]
	yy.A <- 1/pcoef[, "Vmax"] * (1 + pcoef[, "KmA"]/A.uniq[, "[A]"])
	xxplus.A <- c(A.int, xx.A)
	yyplus.A <- c(0, yy.A)
	xxplus.B <- c(B.int, xx.B)
	yyplus.B <- c(0, yy.B)
	xran <- range(c(B.int * 5, A.int, xx.B * 5, xx.A))
	yran <- range(c(0, yy.B, yy.A))
	plot(0, 0, xlim = xran, ylim = yran, las = 1, type = "n", xlab = "", 
		ylab = "")
	axis(2, labels = F, tick = T, pos = 0)
	axis(1, labels = F, tick = T, pos = 0)
	axis(2, labels = F, tick = T, tck = 1, at = 0)
	axis(1, labels = F, tick = T, tck = 1, at = 0)
	title(las = 0, ylab = "1/Vapp", xlab = "(o , 1/[A]); (+ , 1/[B] x 5)")
	#axes(axes = F, las = 0, ylab = "1/Vapp")
	#axes(axes = F, las = 0, xlab = "(o , 1/[A]); (+ , 1/[B] x 5)")
	points(xx.B * 5, yy.B, pch = "+")
	lines(xxplus.B * 5, yyplus.B)
	points(xx.A, yy.A, pch = "o")
	lines(xxplus.A, yyplus.A)
	invisible()
}