function (y=SIR,n=c(6,6,6),Nsim=10000){#try Nsim=10000 first for speed
# This code assumes that SIR is a vector of length 18 containing the SIR 
# results for Flux X, Flux Y, Flux Z, in that order.
F.obs=n[1]*mean(y[1:n[1]])^2+n[2]*mean(y[n[1]+
    1:n[2]])^2+n[3]*mean(y[n[1]+n[2]+1:n[3]])^2
N=18
tr=3
Ybar=mean(SIR)
SST=sum((SIR-Ybar)^2)
Fstat.obs=((N-tr)/(tr-1))*(F.obs-N*Ybar^2)/(SST-F.obs+N*Ybar^2)
F.eq=rep(0,Nsim)
for(i in 1:Nsim){
ind=sample(1:N)
F.eq[i]=n[1]*mean(y[ind[1:n[1]]])^2+
    n[2]*mean(y[ind[n[1]+1:n[2]]])^2+n[3]*mean(y[ind[n[1]+n[2]+1:n[3]]])^2
}
out=hist(F.eq,nclass=100,main="Simulated Randomization Distribution",
   xlab="F-equivalent Test Statistic",col=c("blue","orange"))
abline(v=F.obs,col="red",lwd=2)
pval=mean(F.eq>=F.obs)
text(F.obs+.2,.27*max(out$counts),
    paste("F-equivalent test statistic ",format(signif(F.obs,5))),adj=0)
text(F.obs+.2,.2*max(out$counts),paste("p-value =",format(signif(pval,4))),adj=0)
text(F.obs+.2,.13*max(out$counts),paste("based on ",Nsim," simulations"),adj=0)
readline("hit return/enter\n")
x=seq(0,10,.001)
y=df(x,tr-1,N-tr)
Fstat=((N-tr)/(tr-1))*(F.eq-N*Ybar^2)/(SST-F.eq+N*Ybar^2)
sel=Fstat<10
out=hist(Fstat[sel],breaks=seq(0,10,.2),main="Simulated Randomization Distribution",
xlab="F-Statistic",col=c("blue","orange"),probability=T,ylim=c(0,max(y)))
abline(v=Fstat.obs,col="red",lwd=2)
pvalF=1-pf(Fstat.obs,tr-1,N-tr)
text(Fstat.obs+.2,.27*max(out$density),
    paste("F-statistic ",format(signif(Fstat.obs,5))),adj=0)
text(Fstat.obs+.2,.2*max(out$density),paste("p-value =",format(signif(pval,4)),
" & p-value =",format(signif(pvalF,4))," from F-distribution"),adj=0)
text(Fstat.obs+.2,.13*max(out$density),paste("based on ",Nsim," simulations"),adj=0)
lines(x,y,col="purple")
legend(0,1,"superimposed F-density",lty=1,col="purple",bty="n")
c(F.obs,pval,Fstat.obs)
}