simpfit <- function(ff, p, 
	pfactor = NULL, ndata = NULL, data = NULL, pnames = NULL, 
	iter = 0, verbose = 1, exittest = 1e-08, prtcycle = 30, quadtest = -1,
	maxquadskip = 4)
{
	if(!is.function(ff))
		stop("simpfit: first argument must be a function")
	if(!is.numeric(p))
		stop(
			"simpfit: second argument must be a numeric vector or matrix"
			)
	if(!is.matrix(p)) {
		p <- matrix(p, nrow = 1, byrow = T, dimnames = list(1, names(
			p)))
	}
	nparm <- ncol(p)
	if(length(pnames) == nparm) {
		dimnames(p) <- list(1:nrow(p), pnames)
	}
	else {
		pnames <- eval(as.expression(dimnames(p))[2])
	}
	if(dim(p)[1] == 1) {
		nparm <- length(p)
		if(length(pfactor) == 0)
			pfactor <- 2
		if(length(pfactor)!=nparm)
			pfactor <- (c(pfactor, rep(pfactor, nparm))[1:nparm])
		pp <- matrix(p, nparm, nparm, byrow = T)
		pp[row(pp) == col(pp)] <- pp[row(pp) == col(pp)] * pfactor
		for(i in 1:nparm) {
			if(pfactor[i]!=1)
				p <- rbind(p, pp[i,  ])
		}
	}
	storage.mode(p) <- "double"
	nvert <- nrow(p)
	dimnames(p) <- list(1:nvert, pnames)
	p.start <- p
	maxiter <- iter
	if(quadtest == -1)
		quadtest <- prtcycle
	qparmndx <- rep(0, nvert - 1)
	hessian <- qmat <- matrix(0, nvert - 1, nvert - 1)
	storage.mode(qmat) <- "double"
	storage.mode(hessian) <- "double"
	stddev <- matrix(0, 1, nvert - 1)
	pcent <- matrix(0, 1, nparm)
	storage.mode(pcent) <- "double"
	storage.mode(stddev) <- "double"
	dimnames(pcent) <- list(NULL, pnames)
	if(length(data)!=0) {
		assign("Udata", data, frame = 1)
		assign("SIMPFIT.temp", data, where = 1)
		if(length(ndata) == 0)
			ndata <- nrow(data)
		f.check <- function(x)
		{
			x <- U(x, Udata)
			if(!is.numeric(x))
				stop("Need a numeric result")
			as.double(x)
		}
	}
	else {
		if(length(ndata) == 0)
			ndata <- nvert
		f.check <- function(x)
		{
			x <- U(x)
			if(!is.numeric(x))
				stop("Need a numeric result")
			as.double(x)
		}
	}
	assign("U", ff, frame = 1)
	z <- .C("simpS",
		Sfunction = list(f.check),
		p = p,
		nparm = as.integer(nparm),
		nvert = as.integer(nvert),
		ndata = as.integer(ndata),
		ndatval = as.integer(0),
		iter = as.integer(iter),
		verbose = as.integer(verbose),
		prtcycle = as.integer(prtcycle),
		maxquadskip = as.integer(maxquadskip),
		qparmndx = as.integer(qparmndx),
		exittest = as.double(exittest),
		quadtest = as.double(quadtest),
		covar = qmat,
		pcent = pcent,
		stddev = stddev,
		pcentval = as.double(0),
		rmsdata = as.double(0),
		hessian = hessian)
	if(length(pnames) == nparm) {
		dimnames(z$stddev) <- list(NULL, pnames[z$qparmndx + 1])
		dimnames(z$covar) <- list(pnames[z$qparmndx + 1], pnames[z$
			qparmndx + 1])
		dimnames(z$hessian) <- dimnames(z$covar)
	}
	cov.unscaled <- z$covar/z$rmsdata^2
	if(!is.null(data)) {
		if(exists("SIMPFIT.temp"))
			data <- SIMPFIT.temp
	}
	ans <- list(model = list(ff), iter = z$iter, coef = z$pcent, std.err = 
		z$stddev, std.dev = z$rmsdata, cov.unscaled = cov.unscaled,
		hessian = z$hessian, start.simplex = p.start, final.simplex = z
		$p, control.values = list(nparm = nparm, nvert = nvert, ndata
		 = ndata, iter = maxiter, exittest = exittest, verbose = 
		verbose, prtcycle = prtcycle, quadtest = quadtest, maxquadskip
		 = maxquadskip, free.param = z$qparmndx), data = data)
	ans
}