function(program = "two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr",
	jump = T)
{
# set device and margin before calling function
#	x11(geometry = "640x460", close = T)
#	par(mar = c(8, 7, 6, 4) + 0.1)
#	mmm.tau <<- numeric()
#
# set up for unix call
	if(jump) {
		program <- paste(program, ".jump", sep = "")
		expt <- "jump"
		EXPT <- "JUMP"
	}
	else {
		program <- paste(program, ".mix", sep = "")
		expt <- "mix"
		EXPT <- "MIXING"
	}
#
# 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])))
	}
	else {
		mmm.tau <<- rbind(mmm.tau, as.numeric(c(mmm.dat[1, "S"], 
			xxx.lab[ii + 3])))
	}
#
# calculate total E and S concentrations, for use in labels
	Etotal <- 0
	Stotal <- 0
	for(i in 2:mcol) {
		if(mmm.col[i]!="S")
			Etotal <- Etotal + mmm.dat[1, i]
		if(mmm.col[i]!="E")
			Stotal <- Stotal + mmm.dat[1, i]
	}
	Factor <- 1
	while(Stotal > (Etotal * 2 * Factor)) Factor <- Factor * 2	#
# set up the main label, bottom label, and axis labels
# note xlabel == botlabel
	mainlabel <- paste(EXPT, " ", mmm.lab[1], "\n", mmm.lab[2], "\n", sep
		 = "")
	botlabel <- paste(mmm.lab[4], "\n", mmm.lab[5], "\n", mmm.lab[6], "\n",
		"[E]total = ", signif(Etotal, 5), "   [S]total = ", signif(
		Stotal, 5), sep = "")
	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, cex = 0.667, 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()
}