"two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr" T function(program = "two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr",
	jump = T)
{
	x11(geometry = "640x460", close = T)
	par(mar = c(7, 6, 4, 4) + 0.1)
	mmm.tau <<- numeric()
	if(jump) {
		program <- paste(program, ".jump", sep = "")
		expt <- "jump"
	}
	else {
		program <- paste(program, ".mix", sep = "")
		expt <- "mix"
	}
#
# run the simulation
	unix(program)	# bring over the results
	mmm.col <- scan(paste("scr.", expt, ".col", sep = ""), what = character(
		))
	mcol <- length(mmm.col)
	mspecies <- mcol - 1
	mmm.dat <- matrix(scan(paste("scr.", expt, ".dat", sep = "")), ncol = 
		mcol, byrow = T)
	mlen <- length(mmm.dat[, 1])	# put column headers into data matrix
	dimnames(mmm.dat) <- list(rep("", mlen), mmm.col)	#
# set up the x and y axis ranges
# note x axis = time = (0,1), with data scaled according to time range
# time ranges == max time plotted, are hardwired in trange, 
# for 4 time scales, i.e., for 4 representations of data on each plot
# tmlen is the number of rows of data, up to the cutoff time in trange
	yrange <- range(mmm.dat[2:mlen, 2:mspecies])
	xrange <- range(mmm.dat[2:mlen, 1])	#
#	trange <- c(0.01, 0.001, 0.0001, 1e-05)
#	tmlen <- numeric()
#	for(i in 1:4) {
#		tmlen <- c(tmlen, length(mmm.dat[mmm.dat[, 1] < trange[i], 
#			1]))
#	}
#
# bring over file with 2 lines of main labels
# and a line flagged EBAR+SBAR, with e+s final values and 1/tau
# then two lines of bottom labels
	mmm.lab <- scan(paste("scr.", expt, ".lab", sep = ""), what = character(
		), sep = "\n")	#
# store [S] vs 1/tau for mixing expts
# store ([E] + [S]) vs 1/tau for jump expts
	xxx.lab <- scan(paste("scr.", expt, ".lab", sep = ""), what = character(
		))
	ii <- match("EBAR+SBAR", xxx.lab)
	if(jump) {
		mmm.tau <<- rbind(mmm.tau, as.numeric(c(xxx.lab[ii + 1], 
			xxx.lab[ii + 3])))
		expt <- "JUMP"
	}
	else {
		mmm.tau <<- rbind(mmm.tau, as.numeric(c(mmm.dat[1, "S"], 
			xxx.lab[ii + 3])))
		expt <- "MIXING"
	}
#
# set up the main label, bottom label, and axis labels
# note xlabel == botlabel
	mainlabel <- paste(expt, mmm.lab[1], "\n", mmm.lab[2])
	botlabel <- paste(mmm.lab[4], "\n", mmm.lab[5], "\n", mmm.lab[6])
	ylabel <- "concentration (M)"
	xlabel <- "log (seconds)"	#
# set up the plot axes and labels, curves and points to be drawn later
	plot(mmm.dat[2:mlen, 1], mmm.dat[2:mlen, 2], log = "x", type = "n",
		pch = ".", las = 1, lab = c(5, 4, 6), xlim = xrange, ylim = 
		yrange, ylab = "", xlab = "")
	axes(axes = F, mgp = c(3, 1, 0), ylab = "", xlab = xlabel)
	axes(axes = F, mgp = c(5, 1, 0), ylab = ylabel, xlab = "")
	title(cex = 0.667, main = mainlabel)
	axis(1, labels = F, tck = 1, lty = 2)	# add a dotted grid
	axis(2, labels = F, tck = 1, lty = 2)	#
# plot the curves
	for(i in 2:mspecies) {
		points(mmm.dat[2:mlen, 1], mmm.dat[2:mlen, i], type = "l")
	}
	invisible()
}
"two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr" T function(program = "two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr",
	jump = T)
{
#	x11(geometry = "640x307", close = T)
#	par(mar = c(8, 7, 4, 4) + 0.1)
#	mmm.tau <<- numeric()
	if(jump) {
		program <- paste(program, ".jump", sep = "")
		expt <- "jump"
	}
	else {
		program <- paste(program, ".mix", sep = "")
		expt <- "mix"
	}
#
# run the simulation
	unix(program)	# bring over the results
	mmm.col <- scan(paste("scr.", expt, ".col", sep = ""), what = character(
		))
	mcol <- length(mmm.col)
	mmm.dat <- matrix(scan(paste("scr.", expt, ".dat", sep = "")), ncol = 
		mcol, byrow = T)
	mlen <- length(mmm.dat[, 1])	# put column headers into data matrix
	dimnames(mmm.dat) <- list(rep("", mlen), mmm.col)	#
# set up the x and y axis ranges
# note x axis = time = (0,1), with data scaled according to time range
# time ranges == max time plotted, are hardwired in trange, 
# for 4 time scales, i.e., for 4 representations of data on each plot
# tmlen is the number of rows of data, up to the cutoff time in trange
	yrange <- range(mmm.dat[, 2])
	xrange <- c(0, 1)
	trange <- c(0.01, 0.001, 0.0001, 1e-05)
	tmlen <- numeric()
	for(i in 1:4) {
		tmlen <- c(tmlen, length(mmm.dat[mmm.dat[, 1] < trange[i],
			1]))
	}
#
# bring over file with 2 lines of main labels
# and a line flagged EBAR+SBAR, with e+s final values and 1/tau
# then two lines of bottom labels
	mmm.lab <- scan(paste("scr.", expt, ".lab", sep = ""), what = character(
		), sep = "\n")	#
# store [S] vs 1/tau for mixing expts
# store ([E] + [S]) vs 1/tau for jump expts
	xxx.lab <- scan(paste("scr.", expt, ".lab", sep = ""), what = character(
		))
	ii <- match("EBAR+SBAR", xxx.lab)
	if(jump) {
		mmm.tau <<- rbind(mmm.tau, as.numeric(c(xxx.lab[ii + 1], 
			xxx.lab[ii + 3])))
		expt <- "JUMP"
	}
	else {
		mmm.tau <<- rbind(mmm.tau, as.numeric(c(mmm.dat[1, "S"], 
			xxx.lab[ii + 3])))
		expt <- "MIXING"
	}
#
# set up the main label, bottom label, and axis labels
# note xlabel == botlabel
	mainlabel <- paste(expt, mmm.lab[1], "\n", mmm.lab[2])
	botlabel <- paste(mmm.lab[4], "\n", mmm.lab[5], "\n", mmm.lab[6])
	ylabel <- "delta [X]"
	xlabel <- paste("s/", trange[1], ", s/", trange[2], ", s/", trange[
		3], ", s/", trange[4])
	xlabel <- botlabel	#
# set up the plot axes and labels, curves and points to be drawn later
	plot(0, 0, type = "n", pch = ".", las = 1, lab = c(5, 4, 6), xlim = 
		xrange, ylim = yrange, ylab = "", xlab = "")
	axes(axes = F, mgp = c(3, 1, 0), ylab = "", xlab = xlabel)
	axes(axes = F, mgp = c(5, 1, 0), ylab = ylabel, xlab = "")
	title(cex = 0.667, main = mainlabel)
	axis(1, labels = F, tck = 1, lty = 2)	# add a dotted grid
	axis(2, labels = F, tck = 1, lty = 2)	#
# plot the curves, with a label showing the time range (time cutoff)
	xx <- 0
	for(i in 1:4) {
		points(mmm.dat[1:tmlen[i], 1]/trange[i], mmm.dat[1:tmlen[i],
			2], type = "l")
		text(mmm.dat[tmlen[i] - xx, 1]/trange[i], mmm.dat[tmlen[i] -
			xx, 2], labels = paste("s/", trange[i]), adj = 1)
		xx <- xx + 2
	}
	invisible()
}
progfile T x11(geometry = "640x920", close = T) function(progarray = progfile, jump = T, device = x11(geometry = "640x920",
	close = T))
{
	device
	par(mfrow = c(2, 1))
	par(mar = c(7, 6, 4, 4) + 0.1)	#
#	rm(mmm.tau)
	mmm.tau <<- numeric()
	for(i in 1:length(progarray))
		simprod(progarray[i], jump = jump)	#
#	synchronize()	#save mmm.tau, even if plot bombs
#	rm(xx, yy, temp, xxx)
	xx <<- mmm.tau[, 1]
	yy <<- 1/mmm.tau[, 2]
	if(jump)
		xlabel <- "[E] + [S]"
	else xlabel <- "[S]"
	plot(xx, yy, pch = 1, lab = c(5, 4, 6), las = 1, ylab = "", xlab = "")
	axes(axes = F, mgp = c(3, 1, 0), ylab = "", xlab = xlabel)
	axes(axes = F, mgp = c(5, 1, 0), ylab = "1/tau", xlab = "")
	temp <- nlsfit(twostep.fit)
	lines(spline(xx, yy - temp$residuals))
	xxx <- signif(temp$coef, 4)
	title(cex = 0.667, main = paste("k-2 =", xxx[1], "  k2 =", xxx[2],
		"  k-1/k1 =", xxx[3]))
	rm(xx, yy)
	invisible()
}
two_step -nooutput -E 1.e-5 -S 1.e-4; prep_sfiles scr two_step -nooutput -E 3.1622776602e-5 -S 3.1622776602e-4; prep_sfiles scr two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr two_step -nooutput -E 3.1622776602e-4 -S 3.1622776602e-3; prep_sfiles scr two_step -nooutput -E 1.e-3 -S 1.e-2; prep_sfiles scr two_step -nooutput -E 3.1622776602e-3 -S 3.1622776602e-2; prep_sfiles scr[1] "two_step -nooutput -E 1.e-5 -S 1.e-4; prep_sfiles scr"                    
[2] "two_step -nooutput -E 3.1622776602e-5 -S 3.1622776602e-4; prep_sfiles scr"
[3] "two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr"                    
[4] "two_step -nooutput -E 3.1622776602e-4 -S 3.1622776602e-3; prep_sfiles scr"
[5] "two_step -nooutput -E 1.e-3 -S 1.e-2; prep_sfiles scr"                    
[6] "two_step -nooutput -E 3.1622776602e-3 -S 3.1622776602e-2; prep_sfiles scr"