function(program = "two_step -nooutput -E 1.e-4 -S 1.e-3; prep_sfiles scr.mix"
	)
{
#	x11(geometry = "640x307", close = T)
#	par(mar = c(8, 7, 4, 4) + 0.1)
#	mmm.tau <<- numeric()
	unix(program)
	mmm.col <<- scan("scr.mix.col", what = character())
	mmm.dat <<- matrix(scan("scr.mix.dat"), ncol = length(mmm.col), byrow
		 = T)
	mlen <<- length(mmm.dat[, 1])
	dimnames(mmm.dat) <<- list(rep("", mlen), mmm.col)
	mmm.lab <<- scan("scr.mix.lab", what = character(), sep = "\n")
	xxx.lab <<- scan("scr.mix.lab", what = character())
	ii <- match("EBAR+SBAR", xxx.lab)
# mixing experiments are analyzed as a funtion of [S], not ([E] + [S])
	mmm.tau <<- rbind(mmm.tau, as.numeric(c(xxx.lab[ii + 1], mmm.dat[1, "S"])))
#	mmm.tau <<- rbind(mmm.tau, as.numeric(c(xxx.lab[ii + 1], xxx.lab[ii + 3])))
	mainlabel <<- paste(mmm.lab[1], "\n", mmm.lab[2])
	botlabel <<- paste(mmm.lab[4], "\n", mmm.lab[5], "\n", mmm.lab[6])
	yrange <<- range(mmm.dat[, 2])	
	#	yrange <<- c(mmm.dat[mlen, 2], mmm.dat[1, 2])
	xrange <<- c(0, 1)
	trange <<- c(0.01, 0.001, 0.0001, 1e-05)	
	#	trange <<- c(1e-05, 0.0001, 0.001, 0.01)
	tmlen <<- numeric()
	for(i in 1:4) {
		tmlen <<- c(tmlen, length(mmm.dat[mmm.dat[, 1] < trange[i],
			1]))
	}
	ylabel <<- "delta [X]"
	xlabel <<- paste("s/", trange[1], ", s/", trange[2], ", s/", trange[
		3], ", s/", trange[4])
	plot(mmm.dat[1:tmlen[1], 1]/trange[1], mmm.dat[1:tmlen[1], 2], type = 
		"l", pch = ".", las = 1, lab = c(5, 5, 5), xlim = xrange, ylim
		 = yrange, ylab = ylabel, xlab = botlabel, main = mainlabel)
	axis(1, labels = F, tck = 1, lty = 2)
	axis(2, labels = F, tck = 1, lty = 2)
	points(mmm.dat[1:tmlen[2], 1]/trange[2], mmm.dat[1:tmlen[2], 2], type
		 = "l")
	points(mmm.dat[1:tmlen[3], 1]/trange[3], mmm.dat[1:tmlen[3], 2], type
		 = "l")
	points(mmm.dat[1:tmlen[4], 1]/trange[4], mmm.dat[1:tmlen[4], 2], type
		 = "l")
	xx <- 0
	for(i in 1:4) {
		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
	}
}