t.test.sst <- function(mu, xbar, s=NULL, n=NULL, se=NULL, conf.level=0.95, alternative="two.sided" ) 
{

# Input checking
if( missing("mu")   ) stop("Fatal error: Hypothesized average (mu) must be specified.\n\n")
if( missing("xbar") ) stop("Fatal error: Sample mean (xbar) must be specified.\n\n")
if( missing("n")    ) stop("Fatal error: Sample size (n) must be specified.\n\n")

samples = max( length(xbar),length(s),length(n),length(se) )
  if(samples==1) method="One Sample t-test"
  if(samples==2) method="Two Sample t-test"
  if(samples>2) {
    if(length(mu)>2)   stop("Fatal error: The t-test is only for one or two populations.\n  Your mu parameter has too many values.")
    if(length(xbar)>2) stop("Fatal error: The t-test is only for one or two populations.\n  Your xbar parameter has too many values.")
    if(length(s)>2)    stop("Fatal error: The t-test is only for one or two populations.\n  Your s parameter has too many values.")
    if(length(n)>2)    stop("Fatal error: The t-test is only for one or two populations.\n  Your n parameter has too many values.")
    if(length(se)>2)   stop("Fatal error: The t-test is only for one or two populations.\n  Your se parameter has too many values.")
  }

if(samples==1) { if( length(mu)>1 ) stop("Fatal error: Only one hypothesized mean is allowed.\n\n") }
if(samples==2) { if( length(mu)>1 ) stop("Fatal error: Only one hypothesized difference in means is allowed.\n\n") }

if( length(xbar)!=samples ) {
  xbar=c(xbar,xbar)
  warning("Warning: The sample means were set equal (only one provided).")
  }
if( length(n)!=samples ) {
  n=c(n,n)
  warning("Warning: The sample sizes were set equal (only one provided).")
  }
if(is.null(se)) {
  if(missing("s")) stop("Fatal error: Sample standard deviation (s) \n  or standard error (se) must be specified.\n\n")
  se = s/sqrt(n)
} else {
  s = se*sqrt(n)
}
if( length(s)!=samples ) {
  s=c(s,s)
  se=c(se,se)
  warning("Warning: The standard errors (and standard deviations) were set equal (only one provided).")
  }



# Calculate the test statistic
if(samples==1) {
  TS = (xbar-mu)/se
  df = n-1
  names(mu)   <- "mean"
} else {
  pse = sqrt(sum(se^2))
  TS = -( diff(xbar) - mu ) / pse
  df = (sum(se^2))^2 / ( se[1]^4/(n[1]-1) + se[2]^4/(n[2]-1) )
  names(mu)   <- "difference in means"
}


# Calculate the p-value
if(alternative=="two.sided") {
  p = 1-pt(abs(TS),df=df)
  p = 2*p
}

if(alternative=="less") {
  p = pt(TS,df=df)
}

if(alternative=="greater") {
  p = 1-pt(TS,df=df)
}


# Calculate the confidence interval
if(samples==1) {
  if(alternative=="two.sided") { lcp=xbar+qt( (1-conf.level)/2, df=df ) * se; ucp=xbar-qt( (1-conf.level)/2, df=df ) * se }
  if(alternative=="less")      { lcp=-Inf; ucp=xbar-qt((1-conf.level),df=df )*se }
  if(alternative=="greater")   { lcp=xbar+qt((1-conf.level),df=df)*se; ucp=Inf }
} else {
  if(alternative=="two.sided") {   
    lcp = -diff(xbar)+qt( (1-conf.level)/2, df=df ) * pse
    ucp = -diff(xbar)-qt( (1-conf.level)/2, df=df ) * pse
  }
  if(alternative=="less")      { lcp=-Inf; ucp=-diff(xbar)-qt((1-conf.level),df=df)*pse }
  if(alternative=="greater")   { lcp=-diff(xbar)+qt((1-conf.level),df=df)*pse; ucp=Inf }
}


# Create results
  names(df)   <- "df"
  names(TS)   <- "t"
  names(xbar) <- paste("  Mean",1:samples)
  names(se)   <- paste("  Std Err",1:samples)
  names(s)    <- paste("  Std Dev",1:samples)
  names(n)    <- paste("  n",1:samples, sep="")
  results <- list(
    statistic   = TS, 
    parameter   = df,
    p.value     = p,
    conf.int    = c(lcp, ucp),
    estimate    = c(xbar,s,se,n),
    null.value  = mu,
    alternative = alternative,
    method      = method,
    data.name   = "Summary Data"
  )
  attr(results$conf.int,"conf.level") <- conf.level
  class(results) <- "htest"

return(results)

}