%This is the code that runs the total power comparison between coded and
%uncoded including decoding power for coded.

clear all

PtRange = logspace(log10(0.04),log10(10),50);
RRange = 1.5e9;%logspace(7,log10(7)+9,45); %R<W !!!
eta1 = 4;
alpha1 = 3;
planck = 1.3807e-23;
temp = 300;
Pe = 1e-12;
W = 3e9;
N0 = planck * temp * W;
freq = 60e9;
lambda1 = 3e8 / freq;



r = 1;
theta1 =0;
boxNum = 200;
numReusePts = 200;


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);
Pd = p(1) .* RRange.^2 + p(2).*RRange + p(3);


Duc = zeros(length(PtRange),length(RRange));
DataUC = zeros(length(PtRange),length(RRange));
DcReuse = zeros(length(PtRange),length(RRange),numReusePts);
DataCReuse = zeros(length(PtRange),length(RRange),numReusePts);
Dc = zeros(length(PtRange),length(RRange));
DataC = zeros(length(PtRange),length(RRange));

%DataMaxRUC = zeros(1,length(PtRange));
%DataMaxRC = zeros(1,length(PtRange));


%d>r and d<r cases together with a mod :) -- worst case at theta offset of 0...
for x = 1:length(PtRange)
    PtUC = PtRange(x);
    x
    for y = 1:length(RRange)
        PtC = PtRange(x) - Pd(y);
        
        R = RRange(y);
        PrUC = min(PtUC, PtUC*lambda1^alpha1 / r^alpha1);
        PrC = min(PtC, PtC*lambda1^alpha1 / r^alpha1);
        
        MaxReuseFac = (W/R)*5;
%        reuseRange = linspace(1,MaxReuseFac,numReusePts);

        
        %%%%%%%%%%%%%%%%%%%%%%uncoded case%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        interfGoalUC =2*PrUC/(qfuncinv(Pe)^2) - N0*R/W;
        totalInterf = 0;
        highBnd = 1e6;
        lowBnd = 1e-6;
        
        while highBnd - lowBnd > 1e-4 * lowBnd
            d = lowBnd + (highBnd-lowBnd)/2;
            
            %
                
                s = r*cos(theta1);
                t = r*sin(theta1);
                
                tEff = mod(t,sqrt(3)*d/2);
                rowCount = floor(t/(sqrt(3)*d/2));
                
                if mod(rowCount,2) == 0
                    sEff = mod(s,d);
                    
                else
                    if mod(s,d)>=d/2
                        sEff = mod(s,d) - d/2;
                    else
                        sEff = mod(s,d) +d/2;
                    end
                end
                
                totalInterfTemp = 0;
                
                for j = -boxNum:boxNum
                    if mod(j,2) == 0
                        for i = -boxNum:boxNum
                            term = min (PtUC, PtUC *lambda1^alpha1 / ((i*d - sEff)^2 + (j*sqrt(3)*d/2 - tEff)^2)^(alpha1/2));
                            totalInterfTemp = totalInterfTemp + term;
                        end
                    else
                        for i = -boxNum:boxNum
                            term = min (PtUC, PtUC *lambda1^alpha1 / ((d/2+i*d - sEff)^2 + (j*sqrt(3)*d/2 - tEff)^2)^(alpha1/2));
                            totalInterfTemp = totalInterfTemp + term;
                        end
                    end
                end
                
                totalInterf = totalInterfTemp - min(PtUC, PtUC *lambda1^alpha1 / r^alpha1);
                
                %
            
            if (totalInterf >interfGoalUC)
                lowBnd = d;
            else
                highBnd = d;
            end
        end
        Duc(x,y) = W/R*2/(sqrt(3)*highBnd^2); %W/R is bandwidth compensation for uncoded.
%        DataUC(x,y) = R*Duc(x,y);
        
        %%%%%%%%%%%%%%%%%%%%%%%%coded case%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         
        reuseRange = 2*W/(R/0.8125);
        if (PtC > 0)
            for reuseInd = 1:length(reuseRange)
                reuseFac = reuseRange(reuseInd);
                % assuming Zhengya's 5.5 db SNR.
                
                interfGoalC = min(PtC, PtC*lambda1^alpha1 / r^alpha1) / (max(10^(0.55),(2^(reuseFac*R/W) - 1) )) - N0/(reuseFac);
                totalInterf = 0;
                highBnd = 1e6;
                lowBnd = 1e-6;
                
                while highBnd - lowBnd > 1e-4*lowBnd
                    d = lowBnd + (highBnd-lowBnd)/2;
                    
                    %
                    
                    s = r*cos(theta1);
                    t = r*sin(theta1);
                    
                    tEff = mod(t,sqrt(3)*d/2);
                    rowCount = floor(t/(sqrt(3)*d/2));
                    
                    if mod(rowCount,2) == 0
                        sEff = mod(s,d);
                        
                    else
                        if mod(s,d)>=d/2
                            sEff = mod(s,d) - d/2;
                        else
                            sEff = mod(s,d) +d/2;
                        end
                    end
                    
                    totalInterfTemp = 0;
                    
                    for j = -boxNum:boxNum
                        if mod(j,2) == 0
                            for i = -boxNum:boxNum
                                term = min (PtC, PtC *lambda1^alpha1 / ((i*d - sEff)^2 + (j*sqrt(3)*d/2 - tEff)^2)^(alpha1/2));
                                totalInterfTemp = totalInterfTemp + term;
                            end
                        else
                            for i = -boxNum:boxNum
                                term = min (PtC, PtC *lambda1^alpha1 / ((d/2+i*d - sEff)^2 + (j*sqrt(3)*d/2 - tEff)^2)^(alpha1/2));
                                totalInterfTemp = totalInterfTemp + term;
                            end
                        end
                    end
                    
                    totalInterf = totalInterfTemp - min(PtC, PtC *lambda1^alpha1 / r^alpha1);
                    
                    %
                    
                    if (totalInterf >interfGoalC)
                        lowBnd = d;
                    else
                        highBnd = d;
                    end
                end
                DcReuse(x,y,reuseInd) = reuseFac * 2/(sqrt(3)*highBnd^2);
                DataCReuse(x,y,reuseInd) = R*DcReuse(x,y,reuseInd);
                
            end
            Dc = max(DcReuse,[],3);

            
            [DataC,maxDataCIndex] = max(DataCReuse,[],3);


            if maxDataCIndex == length(reuseRange);
                if abs(DataCReuse(x,y,length(reuseRange)) - DataCReuse(x,y,1)) >=1e-6
                    disp(strcat('R=',num2str(R),', Ptotal=',num2str(PtUC)))
                end
            end

            
            
        end
   
        
        
        
        
    end
    DucMax(x) = max(Duc(x,:));
    DcMax(x) = max(Dc(x,:));
    
    
    
end

% figure, surf(log10(RRange),log10(PtRange),log10(Duc))
% zlabel('Log Density, Uncoded')
% ylabel('Total Power (Pt)')
% xlabel('Rate')
% figure, surf(log10(RRange),log10(PtRange),log10(DataUC))
% zlabel('Log Data Density, Uncoded')
% ylabel('Total Power (Pt)')
% xlabel('Rate')
% 
% figure, surf(log10(RRange),log10(PtRange),log10(Dc))
% zlabel('Log Density, Coded')
% ylabel('Total Power (Pt+Pd)')
% xlabel('Rate')
% figure, surf(log10(RRange),log10(PtRange),log10(DataC))
% zlabel('Log Data Density, Coded')
% ylabel('Total Power (Pt+Pd)')
% xlabel('Rate')
% 
% % for reuseInd = 1:length(reuseRange)
% % figure, surf(log10(RRange),log10(PtRange),log10(DataCReuse(:,:,reuseInd)))
% % zlabel(strcat('Data Density, coded, reuse factor =',num2str(reuseRange(reuseInd))));
% % end
% 
% figure, plot(log10(PtRange),log10(DataMaxRUC))
% hold on, plot(log10(PtRange),log10(DataMaxRC),'r')
% legend('Uncoded','Coded','fontsize',15,'fontname','Times')
% ylabel('Log Max Data Density','fontsize',15,'fontname','Times')
% xlabel('Total Power (Pt+Pd where applicable','fontsize',15,'fontname','Times')
% 
% 
 save dataRuns/triTotalPowerWithPowerGap