function p = PowerDiscreteImperfectCSI(ghat,sigma2w,sigma2e,R,Eta,BER)
gbar = sigma2e/sigma2w;
c1 = 1.6*gbar./(2.^R-1);
c2 = ghat/gbar;
c3 = 0.2*exp(-c2);
p = 1./c1*(c2./(mylambertw(c2./c3*BER))-1);
% y= c3./(c1.*p+1).*exp(c2./(c1.*p+1))

f = @(x) 1/gbar.*exp(-(ghat+x)/gbar).*besseli(0,2/gbar*sqrt(ghat*x));
y = quad(f,0,1);
y2 = 1 - marcumq(sqrt(2*ghat/gbar),sqrt(2*y/gbar));
y - y2