function(k, n0, pstar, conf, rep)
{
# last update 1/30/97
# function to estimate critical values for comparisons with a standard
# when mu0 is known, variances are unknown and unequal.
# returns estimate of (g,h) with conf*100% upper confidence bound on g
#
# k = number of treatments excluding the standard
# n0 = first-stage sample size
# pstar = 1-alpha value (PCS); single value only in this version
# rep = number of replications to use for estimate
# conf = confidence level on upper bound
#
# first get h numerically
	Hstar <- qt(pstar^(1/k), n0 - 1)	#
# then use h to get g
	set.seed(13)
	Z <- matrix(rnorm((k - 1) * rep), ncol = k - 1)
	Zk <- matrix(rnorm(rep), ncol = 1)
	Y <- matrix(rchisq((k - 1) * rep, n0 - 1), ncol = k - 1)
	C <- matrix(rchisq(rep, n0 - 1), ncol = 1)
	Cmat <- matrix(C, ncol = k - 1, nrow = rep, byrow = F)
	denom <- sqrt((n0 - 1) * (1/Y + 1/Cmat))
	G <- sort(apply(cbind(Zk/sqrt(C/(n0 - 1)) + Hstar, Z * denom), 1, max))
	Gstar <- quantile(G, pstar)
	upper <- ceiling(pstar * rep + qnorm(conf) * sqrt(pstar * (1 - pstar) * 
		rep) + 0.5)
	Gupper <- G[upper]
	result <- rbind(Hstar, Gstar, Gupper)
	result
}