% This code takes triangular lattice and shows the density as a function of
% probability of error assuming a saturated transmit power.


clear all

Pt = 1e5;%logspace(-17,3,10);
RRange = 1.5e9;%logspace(7,log10(7)+9,5); %R<W !!!
PeRange = [logspace(-32,-6,10),logspace(-5,-0.5,20)];

alpha1 = 3;
planck = 1.3807e-23;
temp = 300;
Pe = 1e-12;
W = 3e9;
N0 = planck * temp * W;
freq = 60e9;
lambda1 = 3e8 / freq;

GapFromCapacitydB = 3;

GapFromCapacity = 10^(-GapFromCapacitydB/10);


r = 1;
theta1 = 0;
boxNum = 50;
numReusePts = 100;

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

DataMaxRUC = zeros(1,length(PeRange));
DataMaxRC = zeros(1,length(PeRange));

%d>r and d<r cases together with a mod :) -- worst case at theta offset of 0...
for x = 1:length(PeRange)
    x
    Pe = PeRange(x);
    for y = 1:length(RRange)
        
        RUC = RRange(y);
        RC = RRange(y) * (1-(-log2(Pe)*Pe-log2(1-Pe)*(1-Pe)));
        
        
        Pr = min(Pt, Pt*lambda1^alpha1 / r^alpha1);
        MaxReuseFac = (W/RC)*50;
        reuseRange = linspace(1,MaxReuseFac,numReusePts);  

        %%%%%%%%%%%%%%%%%%%%%%uncoded case%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         
%         interfGoalUC =2*Pr/(qfuncinv(Pe)^2) - N0*RUC/W;
%         totalInterf = 0;
%         highBnd = 1e12;
%         lowBnd = 1e-12;
%         
%         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 (Pt, Pt *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 (Pt, Pt *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(Pt, Pt *lambda1^alpha1 / r^alpha1);
%             
%             %
%             
%             if (totalInterf >interfGoalUC)
%                 lowBnd = d;
%             else
%                 highBnd = d;
%             end
%         end
%         Duc(x,y) = W/RUC*2/(sqrt(3)*highBnd^2); %W/R is bandwidth compensation for uncoded.
%         DataUC(x,y) = RUC*Duc(x,y);
%         
%         %%%%%%%%%%%%%%%%%%%%%%%%coded case%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        for reuseInd = 1:length(reuseRange)
            reuseFac = reuseRange(reuseInd);
            
            interfGoalC = GapFromCapacity*min(Pt, Pt*lambda1^alpha1 / r^alpha1) / (2^(reuseFac*RC/W) - 1) - N0/(reuseFac);
            totalInterf = 0;
            highBnd = 1e12;
            lowBnd = 1e-12;
            
            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 (Pt, Pt *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 (Pt, Pt *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(Pt, Pt *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) = RUC*DcReuse(x,y,reuseInd);% the rate corresponds to the actual rate
            % not the rate as seen by the channel coder. That iss why RUC is
            % used.
            
        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(RUC),', Pt=',num2str(Pt)))
            end
        end
       
    end
    DataMaxRUC(x) = max(DataUC(x,:));
    DataMaxRC(x) = max(DataC(x,:));
    DucMaxRUC(x) = max(Duc(x,:));
    DcMaxRC(x) = max(Dc(x,:));
end

% figure, surf(log10(RRange),log10(PeRange),log10(Duc))
% zlabel('Density, Uncoded')
% ylabel('Pe')
% xlabel('Rate')


% 
% figure, surf(log10(RRange),log10(PeRange),log10(DataUC))
% zlabel('Data Density, Uncoded')
% ylabel('Pe')
% xlabel('Rate')
% 
% figure, surf(log10(RRange),log10(PeRange),log10(Dc))
% zlabel('Density, Coded, Pt only','fontsize',16,'fontname','Times')
% ylabel('Pe','fontsize',16,'fontname','Times')
% xlabel('Rate')
% figure, surf(log10(RRange),log10(PeRange),log10(DataC))
% zlabel('Data Density, Coded, Pt only','fontsize',16,'fontname','Times')
% ylabel('Pe','fontsize',16,'fontname','Times')
% xlabel('Rate','fontsize',16,'fontname','Times')
% 
% 



% 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(DucMaxRUC,log10(PeRange))
hold on, plot(DcMaxRC,log10(PeRange),'r')
legend('Uncoded','Coded')
xlabel('Maximum achievable density (nodes/m^2)','fontsize',16,'fontname','Times')
ylabel('log_{10}(Pe)','fontsize',16,'fontname','Times')



save dataRuns/triProbError
