function [] =awgnwaterslide()

clear all;

global Rdata Rdec Rch W boltzmann temperature thermalNoise pathlossexp r lambda1;
global pathloss Enode degree pe;

W = 3e9;
pathlossexp = 3;
r = 10;
lambda1 = 1/200; % 3e8/60e9
boltzmann = 1.38065e-23;
temperature = 300;
Rdata = 1.5e9;
Rdec = Rdata;
Enode = 10e-12;
degree = 4;

Rch = Rdata/W;
thermalNoise = boltzmann * temperature * W;
pathloss = (r/lambda1)^pathlossexp;

log10pevec = [-0.6:-0.25:-5 -5.5:-0.5:-64];%[-0.5:-0.05:-4,-4.25:-0.25:-6,-6.5:-0.5:-10,-11:-1:-40];

for i=1:length(log10pevec)
i

pevec(i) = 10^log10pevec(i);

pe = pevec(i);

term = 1 - hb(pe);

PrStart = thermalNoise *  (     2^(Rch*term)  - 1   );

PShannon(i) = PrStart*pathloss;

[Propt(i),Ptotalmin(i)] = fminbnd(@returnTotalPower,PrStart,PrStart*100,optimset('TolX',1e-40));

%Propt
%Ptotalmin
Ptopt(i) = pathloss * Propt(i) ; 


end
figure;
semilogy(Ptotalmin,log10pevec);
hold
plot(Ptopt,log10pevec,'r')
plot(PShannon,log10pevec,'k')
save dataRuns/awgnwaterslideRun2

%%%%%%%%%%%%%%%%%
function [totalPower] = returnTotalPower(Pr)

global pathloss Rdec Enode Rch degree PrCurrent;

PrCurrent = Pr;

Pt = Pr * pathloss; 

sigmagstart = Pr/ (2^Rch -1);

[sigmagopt, negnmax] = fminbnd(@returnNegnmax,sigmagstart,20*sigmagstart, optimset('TolX',1e-20)); 

nmax = max(-negnmax,1);

iterations = log2((degree-2)*nmax/degree   + 2/degree)   /   log2(degree-1);

Pdec = Enode * Rdec * iterations / Rch;


if Pdec<0
	disp('Pdec<0');
	Pdec
end

totalPower = Pt + Pdec;


%%%%%%%%%%%%%%%%
function [negn] = returnNegnmax(sigmagsq)

global thermalNoise hbinvdeltag pe Rch PrCurrent;

Cg = log2(1+ PrCurrent/sigmagsq);

deltag = 1- Cg/Rch;

hbinvdeltag = hbinv(deltag);

a = 0.5 * (  sigmagsq/thermalNoise - 1 - log( sigmagsq/thermalNoise   )   );

b = ( 1.5 + 2 * log(2/hbinvdeltag)  )  *  (sigmagsq/thermalNoise -1);

c =  log(2* pe/ hbinvdeltag);

n = (-b + sqrt(b^2 - 4 * a * c)  )/ (2 * a);

negn = -n;
