clear all
global log10pe;
global Pt Wused Ractual PtotalRange Rch Rdec;
global RChUse alpha1 planck temp Pe W NO freq lambda1 r theta1 dist p Interference Enode;
alpha1 = 3;
planck = 1.3807e-23;
temp = 300;
Pe = 1e-12;
W = 3e9;
N0 = planck * temp * W;
freq = 60e9;
lambda1 = 3e8 / freq;
Ractual =1.5e9;
Rdec = Ractual;
Enode = 10e-12;

r = 1;
theta1 =0;

PdActual = [144  410 570 770 990]*1e-3;
ThActual = [6.67  16.6 20.8 24.9 29.6]*1e9;
p = polyfit(ThActual,PdActual,2);

p(3) = 0; % forcing the static power to zero.


NoOfBandsRange = linspace(1,6,30);%[1 2 3 4 5 6];% 4 5 6 7 8 10:5:100 150:50:500];%:10:1000];

log10pe =    -12;%[-0.35:-0.05:-2 , -1.99:-0.2:-6, -6:-0.5:-10, -10:-0.5:-32, -32:-1:-64];

DensityRange=[0.02:0.005:0.40  0.405:0.001:.4169];%[0.01:0.03:0.25 0.27:0.02:0.35 0.36:0.001:0.4109];%:0.05:0.25];



for j = 1:length(NoOfBandsRange)
        tic    
    j
    
    NoOfBands = NoOfBandsRange(j);

    InBandDensityRange = DensityRange/NoOfBands;

    Wused = W/NoOfBands;

    
    RChUse = Ractual/Wused;

    for i = 1:length(InBandDensityRange)
    
        InBandDensity = InBandDensityRange(i);
    
        dist = sqrt(2/(sqrt(3)*InBandDensity));
    
        [InterfOptimal(i,j),Ptopt(i,j), Ptotalmin(i,j)] = findptopt();
        
        Optimal_transmit_power = Ptopt(i,j);
        Total_power = Ptotalmin(i,j);
        
        
   
     end
       toc

end
save dataRuns/DensityUpperBound.mat

[OptimalTotalPower,indexes] = min(Ptotalmin,[],2);
figure;
semilogx(OptimalTotalPower,DensityRange);