simulH1<-function(n=5,es=1,sd=1,unil=FALSE,plot=TRUE,power=TRUE,wilcox=FALSE) {
  if(n<2) {
    n<-2;
    cat("Nombre de données a été fixé à 2\n");
  }

  if(es<0) {
    es<-abs(es);
    cat("Différence entre H0 et H1 ne peut pas être négative, valeur absolue prise\n");
  }
  df<-(n-1)*2;
  b<-rnorm(n,0,sd);
  a<-rnorm(n,es,sd);
  my.sd<-sd;
  ncp<-es*(sqrt(n)/my.sd);
  if(plot) {
    par(mfrow=c(2,1),mar=c(5,4,1,1)+0.1);

    vec<-seq(-10,ncp+10,0.001);
    density.H0<-dt(vec[vec<10],df);
    if(es>0) {
      lim.H1<-ncp-4*my.sd;
      density.H1<-dt(vec[vec>=lim.H1],df,ncp);
      my.lim<-c(
        min(vec[vec<10][min(which(density.H0>0.001))],vec[vec>=lim.H1][min(which(density.H1>0.001))]),
        max(vec[vec<10][max(which(density.H0>0.001))],vec[vec>=lim.H1][max(which(density.H1>0.001))])
      );
    }
    else {
      my.lim<-vec[range(which(density.H0>0.001))];
    }
    plot(vec[vec<10],density.H0,type="l",col="red",xlab="Statistique de T",ylab="Densité distribution de T", xlim=my.lim);
    if(es>0) {
      lines(vec[vec>=lim.H1],density.H1,type="l",col="blue");
    }
    if(unil) {
      perc<-qt(0.95,df);
      if(power) {
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
      tt<-t.test(a,b,var.equal=TRUE,alternative="greater");
      wt<-wilcox.test(a,b,alternative="greater");
      vec.1<-seq(perc,my.lim[2],0.001);
      if(es>0) {
        lines(vec.1,dt(vec.1,df,ncp),type="h",col="blue");
      }
      lines(vec.1,dt(vec.1,df),type="h",col="red");
    }
    else {
      tt<-t.test(a,b,var.equal=TRUE);
      wt<-wilcox.test(a,b);
      perc<-qt(0.975,df);
      if(power) {
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
      vec.1<-seq(perc,my.lim[2],0.001);
      if(es>0) {
        lines(vec.1,dt(vec.1,df,ncp),type="h",col="blue");
      }
      lines(vec.1,dt(vec.1,df),type="h",col="red");
      perc<-qt(0.025,df);
      vec.1<-seq(my.lim[1],perc,0.001);
      lines(vec.1,dt(vec.1,df),type="h",col="red");
    }
    pt<-tt$p.value;
    t<-tt$statistic;
    pw<-wt$p.value;
    w<-wt$statistic;
    abline(h=0);
    abline(v=t,col="green");
    if(es>0) {
      height<-max(density.H1)/2;
      x<-(my.lim[2]-my.lim[1])*0.3/13; #the ideal size of the white box behind H1, changing with xlim
      rect(ncp-x,height-0.02,ncp+x,height+0.02,col="white",border="white");
      rect(-x,height-0.02,x,height+0.02,col="white",border="white");
      text(ncp,height,labels="H1",col="blue");
      text(0,height,labels="H0",col="red");
    }
    else {
      height<-max(density.H0)/2;
      x<-(my.lim[2]-my.lim[1])*0.3/13; #the ideal size of the white box behind H1, changing with xlim
      rect(-x,height-0.02,x,height+0.02,col="white",border="white");
      text(0,height,labels="H0",col="red");
    }
    legend(my.lim[1],max(density.H0),legend="T obs.",pch="-",col="green");

#     my.range<-range(qnorm(c(0.0001,0.9999),es,sd/sqrt(n)),tt$conf.int[1]-0.5,tt$conf.int[2]+0.5);
#    vec.1<-seq(my.range[1],my.range[2],0.001);
#    density.mean<-dnorm(vec.1,es,sd/sqrt(n));
#    plot(vec.1,density.mean,type="l",xlab="Distribution de la différence de moyenne",ylab="Densité de probabilité");
#    abline(v=mean(a)-mean(b),col="purple",lty=2);
#    height<-max(density.mean)/2;
#    segments(tt$conf.int[1],height,tt$conf.int[2],height,col="purple");
#    segments(tt$conf.int[1],height+(height*0.1),tt$conf.int[1],height-(height*0.1),col="purple");
#    segments(tt$conf.int[2],height+(height*0.1),tt$conf.int[2],height-(height*0.1),col="purple");

    tta<-t.test(a);
    ttb<-t.test(b);
    my.range<-range(qnorm(0.9999,es,my.sd/sqrt(n)),qnorm(0.0001,0,my.sd/sqrt(n)),tta$conf.int[1]-0.5,ttb$conf.int[1]-0.5,tta$conf.int[2]+0.5,ttb$conf.int[2]+0.5);
    vec.1<-seq(my.range[1],my.range[2],0.001);
    density.a<-dnorm(vec.1,es,my.sd/sqrt(n));
    density.b<-dnorm(vec.1,0,my.sd/sqrt(n));
    plot(vec.1,density.a,type="l",xlab="Distribution de la moyenne de chaque groupe",ylab="Densité de probabilité");
    lines(vec.1,density.b);
    abline(v=mean(a),col="orange",lty=2);
    abline(v=mean(b),col="cyan",lty=2);
    height<-max(density.a)/2;
    segments(tta$conf.int[1],height,tta$conf.int[2],height,col="orange");
    segments(tta$conf.int[1],height+(height*0.1),tta$conf.int[1],height-(height*0.1),col="orange");
    segments(tta$conf.int[2],height+(height*0.1),tta$conf.int[2],height-(height*0.1),col="orange");
    segments(ttb$conf.int[1],height,ttb$conf.int[2],height,col="cyan");
    segments(ttb$conf.int[1],height+(height*0.1),ttb$conf.int[1],height-(height*0.1),col="cyan");
    segments(ttb$conf.int[2],height+(height*0.1),ttb$conf.int[2],height-(height*0.1),col="cyan");
  }
  else {
    if(unil) {
      tt<-t.test(a,b,var.equal=TRUE,alternative="greater");
      wt<-wilcox.test(a,b,alternative="greater");
      if(power) {
        perc<-qt(0.95,df);
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
    }
    else {
      tt<-t.test(a,b,var.equal=TRUE);
      wt<-wilcox.test(a,b);
      if(power) {
        perc<-qt(0.975,df);
        power.t<-1-pt(perc,df,ncp);
      }
      else {
        power.t<-1.0;
      }
    }
    pt<-tt$p.value;
    t<-tt$statistic;
    pw<-wt$p.value;
    w<-wt$statistic;
  }
  if(power) {
    power.w<-power.wilcox.test(es=es,n=n,sd=sd);
  }
  else {
    power.w<-1.0;
  }

  if(wilcox) {
    res<-matrix(c(t,pt,power.t,w,pw,power.w),ncol=3,nrow=2,byrow=TRUE);
    rownames(res)<-c("Test de T","Test de Wilcoxon");
    colnames(res)<-c("Statistiques","P-value","Puissance");
  }
  else {
    res<-matrix(c(t,pt,power.t),ncol=3,nrow=1,byrow=TRUE);
    rownames(res)<-c("Test de T");
    colnames(res)<-c("Statistiques","P-value","Puissance");
  }

  return(res);
}

power.wilcox.test<-function(es=1,n=5,sd=1,unil=FALSE) {
  wilcox.pval<-function(es=1,n=5,sd=1,unil=FALSE) {
    a<-rnorm(n,0,sd);
    b<-rnorm(n,es,sd);
    if(unil) {
      wilcox.test(a,b,alternative="greater")$p.value;
    }
    else {
      wilcox.test(a,b)$p.value;
    }
  }

  pvals<-replicate(1000,wilcox.pval(es,n,sd));
  return(sum(pvals<0.05)/1000);
}