function [load, throughput,optimum_load, optimum_throughput] = intermittentNodesPerf(P_EH,P_Rx,E_Tx,E_EH_mean,E_EH_sigma)
%MAX_THROUGHPUT_LOAD
%% Calculation of throughput versus Load
% Each time a harvesting event happens a node can choose whether
% to transmit or to use the energy for listening instead.
% As the number of packets a node processes increases the
% number of transmissions required increases.
% Consider where all nodes have the same number to transmit
% as recieve. By always transmitting after a harvesting event
% there is very little energy (and therefor time) to receive,
% so no communication can happen.
% Conversely, by always receiving no transmissions occur so there
% is no communication happening.
% There must be a optimum point in between.
%
% The tx load, $R$, is the product of the transmission
% probability, $p(Tx)$, and the energy harvesting rate, $\lambda_{EH}$.
% Assume a harvester, with power $P_{EH}$, is chosen to deliver
% bursts of energy close to that required for transmission,
% $E_{EH}$ at the harvesting rate.
%
% The equation governing the achieveable throughput, $\theta$ at tx load is
% 
% $$\theta = \frac1P_{Rx}\left(P_{EH}R - E_{TX}R^2\right)$$
%
% The load limit is calculated from the probability it harvests over the
% $E_{TX}$
% threshold in any given event.
    % 2% error due to ignoring events below zero at $dev = 0.5$
    p_no_tx = normcdf(1 * E_Tx,E_EH_mean, E_EH_sigma);
    p_lt2_tx = normcdf(2 * E_Tx,E_EH_mean, E_EH_sigma);
    p_1_tx = p_lt2_tx-p_no_tx;
    p_lt3_tx = normcdf(3 * E_Tx,E_EH_mean, E_EH_sigma);
    p_2_tx = p_lt3_tx - p_lt2_tx;
    % At E_{EH}/E_{Tx} = 2 with a deviation ratio = 0.5
    % this actually ignores 2% which are 4*E_{Tx}. Be careful with high
    % deviation ratios or high E_EH_mean
    p_lt4_tx = normcdf(4 * E_Tx,E_EH_mean, E_EH_sigma);
    p_3_tx = p_lt4_tx - p_lt3_tx;
    rate_EH = P_EH/E_EH_mean;
    load_lim = rate_EH * (p_1_tx + p_2_tx*2+p_3_tx*3);
    load = linspace(0,load_lim,100);
    throughput = 1/P_Rx*(P_EH*load-E_Tx*load.^2);
    optimum_load = min(P_EH/2/E_Tx,load_lim);
    optimum_throughput = 1/P_Rx*(P_EH*optimum_load-E_Tx*optimum_load.^2);
end