function (Nsim = 10000) 
{
    par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)
    SIR = flux$SIR
    FLUX = flux$FLUX
    m.SIR = mean(SIR)
    SIR.X = mean(SIR[FLUX == "X"])
    SIR.Y = mean(SIR[FLUX == "Y"])
    mm = min(SIR[FLUX == "X"], SIR[FLUX == "Y"])
    MM = max(SIR[FLUX == "X"], SIR[FLUX == "Y"])
    qqplot(SIR[FLUX == "X"], SIR[FLUX == "Y"], xlab = expression("SIR with Flux X   (" ~ 
        log[10](Ohm) ~ ")"), ylab = expression("SIR with Flux Y   ( " ~ 
        log[10](Ohm) ~ ")"), xlim = c(mm, MM), ylim = c(mm, MM), 
        pch = 16, col = "red")
    abline(0, 1)
    readline("Hit return\n")
    qqplot(SIR[FLUX == "X"], SIR[FLUX == "Y"], xlab = expression("SIR with Flux X   (" ~ 
        log[10](Ohm) ~ ")"), ylab = expression("SIR with Flux Y   ( " ~ 
        log[10](Ohm) ~ ")"), xlim = c(0, 20), ylim = c(0, 20), 
        pch = 16, col = "red")
    abline(0, 1)
    readline("Hit return\n")
    boxplot(SIR ~ FLUX, data = flux, ylab = expression("SIR ( " ~ 
        log[10](Ohm) ~ ")"), xlab = "Flux")
    text(1.5, 12.5, substitute(bar(FLUX)[Y] - bar(FLUX)[X] == 
        flux.diff, list(flux.diff = format(signif(SIR.Y - SIR.X, 
        4)))), adj = 0)
    points(rep(0.95, 9), flux$SIR[flux$FLUX == "X"], col = "red", 
        pch = 16)
    points(rep(1.95, 9), flux$SIR[flux$FLUX == "Y"], col = "blue", 
        pch = 16)
    out.AD = AD.test(flux$SIR[flux$FLUX == "Y"], flux$SIR[flux$FLUX == 
        "X"])
    DIFF = SIR.Y - SIR.X
    readline("Hit return\n")
    randomization.ref.dist = combn(1:18, 9, fun = mean.fun, y = SIR)
    D = 2 * (randomization.ref.dist - m.SIR)
    p1 = mean(D < DIFF)
    p2 = mean(D > -DIFF)
    hist(D, nclass = 200, probability = T, col = c("red", "blue"), 
        main = expression("Randomization Reference Distribution of" ~ 
            bar(SIR)[Y] - bar(SIR)[X]), xlab = expression(bar(Y) - 
            bar(X) == bar(SIR)[Y] - bar(SIR)[X]))
    abline(v = c(DIFF, -DIFF), lwd = 2, col = "orange")
    text(1.05 * DIFF, 0.3, substitute(P(bar(Y) - bar(X) <= yxdiff) == 
        p1, list(yxdiff = signif(DIFF, 4), p1 = signif(p1, 4))), 
        adj = 0.5, srt = 90)
    text(-1.05 * DIFF, 0.3, substitute(P(bar(Y) - bar(X) >= yxdiff) == 
        p2, list(yxdiff = signif(-DIFF, 4), p2 = signif(p2, 4))), 
        adj = 0.5, srt = -90)
    readline("Hit return\n")
    qqnorm(D, ylab = expression(bar(Y) - bar(X)), pch = ".")
    qqline(D, col = "blue")
    readline("Hit return\n")
    W = D/sqrt(sum((SIR - m.SIR)^2))
    t.SIR = sort((W/sqrt(1 - 9 * 9 * W^2/18)) * sqrt(18 - 2)/sqrt(1/9 + 
        1/9))
    N.perm = length(t.SIR)
    t.quant = qt((1:N.perm)/(N.perm + 1), 16)
    plot(t.quant, t.SIR, pch = ".", xlab = expression(t[16] ~ 
        "quantiles"), ylab = "ordered randomization t-statitics")
    abline(0, 1, col = "blue")
    readline("Hit return\n")
    t.sample = sort(rt(N.perm, 16))
    mt = t.sample[1]
    Mt = t.sample[N.perm]
    plot(t.quant, t.sample, pch = ".", xlab = expression(t[16] ~ 
        "quantiles"), ylab = "simulated t-statitics", xlim = c(mt, 
        Mt), ylim = c(mt, Mt))
    abline(0, 1, col = "blue")
    readline("Hit return\n")
    sim.ref.dist = NULL
    for (i in 1:Nsim) {
        SIR.star = sample(SIR)
        sim.ref.dist = c(sim.ref.dist, mean(SIR.star[1:9]))
    }
    D.star = 2 * (sim.ref.dist - m.SIR)
    m = min(D.star, D)
    M = max(D.star, D)
    qqplot(D, D.star, xlab = expression(bar(y) - bar(x) ~ "for all combinations"), 
        , xlim = c(m, M), ylim = c(m, M), ylab = substitute(bar(y) - 
            bar(x) ~ "for all" ~ Nsim ~ "sampled combinations", 
            list(Nsim = Nsim)), pch = ".")
    abline(0, 1, col = "blue")
    p1.hat = mean(D.star < DIFF)
    p2.hat = mean(D.star > -DIFF)
    p1.hat = signif(p1.hat, 4)
    p2.hat = signif(p2.hat, 4)
    p1 = signif(p1, 4)
    p2 = signif(p2, 4)
    text(m, M, substitute(hat(p)[1] == p1.hat, list(p1.hat = p1.hat)), 
        adj = 0, cex = 1)
    text(m, 0.87 * M, substitute(hat(p)[2] == p2.hat, list(p2.hat = p2.hat)), 
        adj = 0, cex = 1)
    text(M, 0.9 * m, substitute(p[1] == p1, list(p1 = p1)), adj = 1, 
        cex = 1)
    text(M, m, substitute(p[2] == p2, list(p2 = p2)), adj = 1, 
        cex = 1)
    out = list(out.AD = out.AD, SIR.X = SIR.X, SIR.Y = SIR.Y, 
        DIFF = DIFF, p1 = p1, p2 = p2)
}