function (nn = c(6, 6, 6), Nsim0 = 10000, Nsim1 = 10000, alpha = 0.05, 
    y = as.vector(as.matrix(Flux3))) 
{
    k = length(nn)
    N = sum(nn)
    XbarYbar = NULL
    for (i in 1:Nsim0) {
        ind = sample(1:N, N)
        XbarYbar = c(XbarYbar, mean(y[ind[1:nn[1]]]) - mean(y[ind[(nn[1] + 
            1):(nn[1] + nn[2])]]))
    }
    tcrit = quantile(abs(XbarYbar), 0.95)
    par(mfrow = c(1, 1))
    hist(XbarYbar, nclass = 101, xlab = expression(bar(X) - bar(Y)))
    Xbar.mat = NULL
    for (i in 1:Nsim1) {
        Xbar.vec = NULL
        ind = sample(1:N, N)
        start = 0
        for (j in 1:k) {
            Xbar.vec = c(Xbar.vec, mean(y[ind[start + (1:nn[j])]]))
            start = start + nn[j]
        }
        Xbar.mat = rbind(Xbar.mat, Xbar.vec)
    }
    rejections = NULL
    for (i in 1:(k - 1)) {
        for (j in (i + 1):k) {
            rejections = cbind(rejections, abs(Xbar.mat[, i] - 
                Xbar.mat[, j]) > tcrit)
        }
    }
    rejections0 = apply(rejections, 1, max)
    rates = apply(rejections, 2, mean)
    rate0 = mean(rejections0)
    c(rates, rate0)
}