function y = BERImperfectCSI(p,ghat,sigma2w,sigma2e,R,BER)
gbar = sigma2e/sigma2w;
% f = @(x) 0.2*exp(-1.6*p*x./(2^R(end)-1))/gbar.*exp(-(ghat+x)/gbar).*besseli(0,2/gbar*sqrt(ghat*x));
% y = quad(f,0,ghat+gbar+10*sqrt(gbar^2+2*ghat*gbar));
c1 = 1.6*gbar./(2.^R-1);
c2 = ghat/gbar;
c3 = 0.2*exp(-c2);
y= c3./(c1*p+1).*exp(c2./(c1*p+1));
p = 1./c1*(c2./(mylambertw(c2./)));