% FOR AWGN
function [InterfOpt, Ptopt, Ptotalmin] = findptopt()

global Interference lambda1 alpha1 r theta1 dist PtotalRange;

boxNum = 100;

d = dist;

PtRange = logspace(-6,1,200);


for PtIndex=1:length(PtRange)

    Pt = PtRange(PtIndex);


    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

    Interference = totalInterfTemp - min(Pt, Pt *lambda1^alpha1 / r^alpha1);

    InterferenceRange(PtIndex) = Interference;
    
    % Interference calculated, now finding total power.

    PtotalRange(PtIndex) = findtotalpowerWithInterference2(Pt);

end


[Ptotalmin,PtotalminIndex] = min(PtotalRange);
Ptopt = PtRange(PtotalminIndex);
InterfOpt = InterferenceRange(PtotalminIndex);
