Simulate power function for a t-test with a contaminated normal distribution

rcn<-function(n,eps,sigmac,u){
ind<-rbinom(n,1,eps)
x<-rnorm(n,u,1)                             #norm with mean=u, sd=1
rcn<-x*(1-ind)+sigmac*x*ind        #produce obs fromcontaminated normal distriubtion with u
return(rcn)                              
}

#  returns a random sample of size n drawn from
#  a contaminated normal distribution with percent
#  contamination eps and variance ratio sigmac


emp_power<-function(nsims,mean){
sigmac<-25
eps<-.25
alpha<-.05
n<-20                               # sample size of 20
tc<-qt(1-alpha,n-1)            #this is 1.729, i.e the reject region,here we suppose u=0
ic<-0
for(i in 1:nsims){
    samp<-rcn(n,eps,sigmac,mean)                    # simulate 20 samples from
                                                                       #the contaminated normal distribution
    ttest<-(sqrt(n)*mean(samp))/var(samp)^.5  # calcuate ts from sample.
    if(ttest > tc){ic<-ic+1}
    }
power<-ic/nsims
return(power)
}

#err<-1.96*sqrt((power*(1-power))/nsims)
#list(empiricalpower=power,error=err)


#We will draw a power function for this simulation
powers<-NULL
u<-seq(0,1,0.001)
for(i in 1:1001){
powers[i]<-emp_power(1000,u[i])
}

plot(u,powers,type="l") # this is power function plot