function Pb=map_bound4(BPH,rho_in_db)
%M=2;
%rho_in_db=7;
%rho=10^(rho_in_db/10);
%k=log(M)/log(2);
%f=exp((-x./(1+x)).*k*rho)./(1-x.^2);
fanout=2^BPH;
u=10^(rho_in_db/10);
u1=BPH*u;
g=inline('exp(namda.*x).*(1/gamma(n)).*x.^(n-1).*exp(-x)','x','n','namda');
h=inline('exp(-namda.*x).*(1/2).*((x./(n.*u)).^((n-2)./4)).*exp((-1/2).*(x+n.*u)).*besseli(n./2-1,sqrt(n.*u.*x))','x','n','u','namda');
for i=0:1:99
    Q(i+1)=quad(g,0,100,[],[],3,i/100);
    Q1(i+1)=quad(g,0,100,[],[],4,i/100);
    J(i+1)=quad(h,0,100,[],[],3,u1,i/100);
    J(i+1)=quad(h,0,100,[],[],4,u1,i/100);
end
P=Q.*J;
P1=Q1.*J1;
y=min((fanout-1).*P+(fanout-1)*(fanout^3-1).*P1);
Pb=(fanout/(2*(fanout-1)))