% How does the expected WU rate affect the optimum goodput
%% Generate throughput load for a node
figure(1);
grid on;
xlabel('Load (s^{-1})')
title({'Variation of optimum load for changing';'number of nodes, n, within wake up range'})
ylabel('Goodput (s^{-1})')
hold on;
P_Rx = 0.35e-6;
P_EH = 0.4e-6;
E_Tx = 0.4e-3;
[load, throughput,optimum_load, optimum_throughput] =...
  intermittentNodesPerf(P_EH, P_Rx, E_Tx, 0.52e-3, 0.2e-3);
  %                   P_EH,   P_Rx,   E_Tx, E_EH_mean, E_EH_sigma
plot(load, throughput);
%% Create a simple model for linking wake ups to transmit rate
n_wu = [1,3,10,33,100];% nodes providing wakeups
T_B = 8e-3;% 1ms of time for each beacon at the receiver
P_wu = 10e-3;%W power consumption of node during wake up event

% Second figure for the wake up prob
figure(1);
for i=1:numel(n_wu)
    p_wu = n_wu(i)*T_B*load;
    adjusted_throughput = throughput./(1+P_wu/P_Rx*p_wu);
    plot(load, adjusted_throughput);
    figure(1);
end
legend_items = num2str([0;n_wu'],'n = %i');
legend(legend_items,'Location', 'NorthWest')

%% With a better model of wakeups 
figure(3)
grid on;
xlabel('Load (s^{-1})')
ylabel('Goodput (s^{-1})')
title({'Variation of optimum load for changing';'density of nodes, \lambda, affecting wake up probability'})
hold on;
T_B = 8e-3;% 1ms of time for each beacon at the receiver
P_wu = 10e-3;%W power consumption of node during wake up event
p_wu = [];

% Using Kouyaza WU Calc
samples_n = 3;
max_lam = 5e-2;
min_lam = 1e-5;
lam = linspace(min_lam,max_lam,samples_n); %lambda
lam = [1e-5,0.01,0.04];

alp = 2.6; %alpha
d = T_B*load; %delta
a = 1; %RF-DC Conversion efficiency
P_db =15; % BS Tx power in dBm
P = 10.^((P_db/10)-3);
T_db =-20; % BS Tx power in dBm
T = 10.^((T_db/10)-3);
for i = 1:numel(d)
    [~,~, p_wu(i,:), ~, ~] = wake_bounds(lam,d(i),alp,a,P,T);
end
% Calc adjusted Throughput
adjusted_throughput = repmat(throughput,samples_n,1)./(1+P_wu/P_Rx*p_wu');
plot(load, adjusted_throughput);
grid on
legend_items = num2str(lam','\\lambda = %.2f');
legend(legend_items,'Location', 'NorthWest')

%% Find Maxima of functions with Optimisation Toolbox
samples_n = 21;
max_lam = 5e-2;
min_lam = 1e-5;
lam = linspace(min_lam,max_lam,samples_n); %lambda
R_lb = 0;
R_ub = optimum_load;
R_max=[];
G_max=[];
G_unadjusted=[];
optimset('TolX',R_ub/1e6);
for i = 1:numel(lam)
    goodput_function = @(R_o)-1/P_Rx*(P_EH*R_o-E_Tx*R_o.^2)/...
        (1+P_wu/P_Rx*p_wu_integral(lam(i),T_B*R_o,alp,a,P,T));
    [R_max(i), neg_G_max] = fminbnd(goodput_function, R_lb,R_ub,optimset('TolX',R_ub/1e6));
    G_unadjusted(i) = -goodput_function(R_ub);
    G_max(i) = -neg_G_max;
end
l=plot(R_max([1 5 13]),G_max([1 5 13]),'o');
set(l, 'DisplayName', 'Max Rate');

figure(6);
grid on;
title('Maximum possible Goodput and the rate to give that value');
xlabel('Density of Nodes, \lambda');
ylabel('Throughput/Goodput (s^{-1})');
hold on;
plot(lam, R_max, '-o');
plot(lam, G_max, '-o');