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(7, 6, 4, 4) + 0.1)
#
# set up for unix call
	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 scale [S]: ([E]total <= [S]total <= 2 * [E]total)
# later, print [S] scale factor on plot
	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
	mmm.scaled <- mmm.dat
	mmm.scaled[, "S"] <- mmm.scaled[, "S"]/Factor
	yrange <- range(mmm.scaled[2:mlen, 2:mcol])
	xrange <- range(mmm.scaled[2:mlen, 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")	#
#
# set up the main label, bottom label, and axis labels
# note xlabel == botlabel
	mainlabel <- paste(expt, ":  [E]total = ", signif(Etotal, 5), 
		"   [S]total = ", signif(Stotal, 5), "\n", mmm.lab[4], sep = ""
		)
	botlabel <- paste("\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.scaled[2:mlen, 1], mmm.scaled[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, sub = botlabel, 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
	xx <- 2
	for(i in 2:mcol) {
		points(mmm.scaled[2:mlen, 1], mmm.scaled[2:mlen, i], type = "l"
			)
		if(mmm.col[i] == "S")
			species.lab <- paste(mmm.col[i], "/", as.character(
				Factor), sep = "")
		else species.lab <- mmm.col[i]
		text(mmm.scaled[mlen - xx, 1], mmm.scaled[mlen - xx, i], labels
			 = species.lab, adj = 1)
		xx <- xx + 10
	}
	invisible()
}
