z.test <- function(x, y=NULL, 
		sigmax, 
		sigmay=sigmax, 
		mu = 0, 
		alternative="two.sided", 
		conf.level=0.95
	) {

# Test for good input

error <- ""

if( sigmax <= 0 || sigmay<=0) { 
	error="Error: Variance out of bounds" 
	}
if(0>conf.level || 1<conf.level) { 
	error="Error: Confidence interval out of bounds"
	}
if(alternative!="two.sided" && alternative!="less" && alternative!="greater") {
	error = "Error: Alternative specification not recognized"
	}

if(error!="") {
return(error)
}



if(is.null(y)) {
  method <- "One-sample z-test"
  sx     <- sqrt( sigmax^2/length(x) )
  xbar   <- mean(x)
  z      <- (xbar-mu)/se
  est    <- xbar
  names(est) <- paste("mean of",deparse(substitute(x)))
  dn     <- deparse(substitute(x))
  nv     <- "mean"
} else {  
  method <- "Two-sample z-test"
  sx     <- sqrt( sigmax^2/length(x) )
  sy     <- sqrt( sigmay^2/length(y) )
  se     <- sqrt( sx^2 + sy^2 )
  xbar   <- mean(x)
  ybar   <- mean(y)
  z      <- ((xbar-ybar)-(mu))/se
  est    <- c(xbar,ybar)
  names(est)[1] <- paste( "mean of",deparse(substitute(x)) )
  names(est)[2] <- paste( "mean of",deparse(substitute(y)) )
  dn     <- paste(deparse(substitute(x)),"(s=",sigmax,") and ",deparse(substitute(y)),"(s=",sigmay,")",sep="")
  nv     <- "difference in means"
  if(sigmax!=sigmay) method <- paste(method,", unequal variances",sep="")
}

e <- sum(est)

# Particular info
alpha  <- 1-conf.level
if(alternative=="two.sided") {
  a   <- pnorm(abs(z)) - pnorm(-abs(z))
  lcp <- e - qnorm(1-alpha/2) * se
  ucp <- e + qnorm(1-alpha/2) * se
}
if(alternative=="less") {
	a   <- 1 - pnorm(z)
	lcp <- -Inf
	ucp <- e + qnorm(1-alpha) * se
}
if(alternative=="greater") {
	a   <- pnorm(z)
	lcp <- e - qnorm(1-alpha) * se
	ucp <- Inf
}


parmtr        <- Inf
names(parmtr) <- "df"
stat          <- z
names(stat)   <- "z"
p             <- 1-a
names(mu)     <- nv

# Create results
results <- list(
	statistic=stat, 
	parameter = parmtr,
	p.value=p,
	conf.int=c(lcp, ucp),
	estimate=est,
	null.value=mu,
	alternative=alternative,
	method=method,
	data.name=dn,
	sigmax=sigmax,
	sigmay=sigmay
	)
attr(results$conf.int,"conf.level") <- conf.level
class(results) <- "htest"

return(results)
}