function y = PowerImperfectCSIApprox(p,ghat,gbar,g0)
% x = linspace(0,ghat+gbar+10*sqrt(gbar^2+2*ghat*gbar),500);
% g = 1./gbar.*exp(-(ghat+x)/gbar).*besseli(0,2/gbar*sqrt(ghat*x));
% y = 2/gbar*ncx2pdf(2*x/gbar,2,2*ghat/gbar);
K = ghat/gbar;
m = (K+1)^2/(2*K+1);
b = m/(ghat+gbar);
% Omegap = sigma2w*(ghat+gbar);
% b = Omegap/m;
% yapp = sigma2w*gampdf(sigma2w*x,m,b);
% plot(x,g,'b',x,y,'r--',x,yapp,'g-.')
% f = @(x) x./(1+p*x)/gbar.*exp(-(ghat+x)/gbar).*besseli(0,2/gbar*sqrt(ghat*x));
% f =  @(x) x./(1+p*x)*2/gbar.*ncx2pdf(2*x/gbar,2,2*ghat/gbar);
% y = quad(f,0,ghat+gbar+10*sqrt(gbar^2+2*ghat*gbar))-lambda/(w/log(2));
% Gamma_a_x = gamma(-m)*(1-gammainc(b/p,-m));
Gamma_a_x = maple('GAMMA',double(-m),double(b/p);
y = (b/p)^(m+1)*exp(b/p)*Gamma_a_x-g0;