PowerSim=function(m=7,n=7,upper=T,alpha=.05,
DeltaSigs=c(.2,.4,.6,.8,1,1.5,2),Nsim=1000){
out=NULL
mx=n*(n+1)/2
N=m+n
p1=pnorm(DeltaSigs/sqrt(2))
if(!exists("pmnorm")) library(mnormt)
p2=NULL
for(i in 1:length(DeltaSigs)){
  zx = DeltaSigs[i]/sqrt(2)
  p2[i]=pmnorm(c(zx,zx),c(0,0),varcov=matrix(c(1,.5,.5,1),ncol=2))
}
EWXY=m*n*p1
VarWXY=m*n*(p1*(1-p1)+(N-2)*(p2-p1^2))
DeltaSig=0
out=combn(1:N,n,FUN=sum)
out=out-mx
W.u= sort(unique(out))
power=NULL
i=0
for(wu in W.u){
  i=i+1
  power[i]=mean(out>=wu)
}
sel0=power<=alpha
sel1=power>alpha
power0=max(power[sel0])
power1=min(power[sel1])
c0=min(W.u[sel0])
c1=max(W.u[sel1])
if(upper==T){
  alphax=power0
  cx=c0}else{
  alphax=power1
  cx=c1
}
wu=c(c0,c1)
powerAlpha=c(power0,power1)
Pow1=1-pnorm((cx-.5-EWXY)/sqrt(VarWXY))
Pow2=pnorm(sqrt(3*m*n/((N+1)*pi))*DeltaSigs-qnorm(1-alphax))
Pow3=pnorm(sqrt(3*m*n/((N)*pi))*DeltaSigs-qnorm(1-alphax))
Pow1.0=1-pnorm((cx-.5-m*n/2)/sqrt(m*n*(N+1)/12))

powerx=NULL
j=0
for(DeltaSig in DeltaSigs){
   j=j+1
   out=NULL
   for(i in 1:Nsim){
      x=rnorm(m,0,1)
      y=rnorm(n,DeltaSig,1)
      rx=rank(c(y,x))
      out[i]=sum(rx[1:n])-mx
   }
   powerx[j]=mean(out>=cx)
}
CxAlpha=round(cbind(c(c0,c1),powerAlpha),3)
Power=rbind(c(0,DeltaSigs),c(alphax,powerx),c(Pow1.0,Pow1),c(alphax,Pow2),c(alphax,Pow3))
Power=round(Power,3)
list(CxAlpha=CxAlpha,Power=Power,powerx=c(alphax,powerx))
}