function [P_Fu,P_Fl, P_W, P_Sl, P_Su] = wake_bounds(lam_arr,d,alp,a,P,T,varargin)
%WAKE_PROB Calculate probability of wake up and bounds for false wake up
%   lam = density of BS, alp = path loss, d = transmission duty cycle,
%   a = Rectenna efficiency, P = Tx Power, T= Wake Threshold Power
    p = inputParser;
    addOptional(p,'Tidy','Raw',@ischar);
    parse(p,varargin{:});
    for i = 1:numel(lam_arr)
        lam = lam_arr(i);
        if(max(d.*lam)>0.24 || min(d.*lam)<1e-5 ...
         ||max(alp) > 4.0 || min(alp) < 2.3...
         ||max(P./a./T)>10^10 || min(P./a./T) <10)
            warning('Out of bounds for accurate integration of wake up prob.');
        end
        Pn_rtalp = (a*P./T).^(1/alp); %m, radius of wake up
        lpPn2a = lam*pi*Pn_rtalp.^2;%Orange
        P_rtalp = (a*P)^(1/alp); %m
        lpP2a = lam*pi*P_rtalp^2;
        F_exp = exp(-lpPn2a);
        P_Fu = 1-F_exp...
                +1/(1-d)*(1-(2*d*lpPn2a)/(alp-2))...
                        .*(F_exp-F_exp.^d);
        P_Fl = 1 - F_exp+ 1/(1-d).*(F_exp-F_exp.^d);
        f_arg = @(u) 2*pi*d/alp*u.^(2/alp)*lpP2a;
        integrand = @(u)1/pi./u.*exp(-u.*T...
        -f_arg(u)/(tan(2*pi/alp)))...
        .*sin(f_arg(u));
    %   P_W(i) = [integral(integrand, 0, 4/(alp-2)^3*(lam*d)^0.6*1e8)];
        P_W(i) = [integral(integrand, 0, 3*(P/a/T)^(0.7)/(alp-2)^3.3*(lam*d)^0.6*1e7)];
        if(strcmp(p.Results.Tidy,'Tidy'))
            if(P_W(i)>1||P_W(i)<-1)
                P_W(i:numel(lam_arr)) = ones([1,numel(lam_arr)-i+1]);
                break;
            end
        end
        
    end
    %% Make limited version of upper and lower limits
    P_Fulim = min(P_Fu, 1);
    P_Fllim = max(P_Fl, 0);
    P_Sl = 1/d*P_W - (1-d)/d*P_Fulim;
    P_Su = 1/d*P_W - (1-d)/d*P_Fllim;
end