function [n] = returnnforsigmag2Density(sigmagsq)  % Remember you have to specify Pt as well, and make sure sigmagsq > Pt/Pshannon

global TotalNoise Preceived;
global log10pe RChUse;
global sigmagsq2

max_iter = 30;
sigmagsq2 = sigmagsq;



Cg = 0.5*log2(  1  +  Preceived/sigmagsq  );

deltag = 1 - Cg/RChUse;

if deltag<0
    
    disp('delta(G) negative');
    deltag
end

logpe = log10pe* log(10);

hbinvdeltag = hbinv(deltag);

%epsilon =  (3 + 4*log(2) + 4*log(1/hbinvdeltag) )*sigmagsq;  

%epsilon fixed now.

%epsilonprime = epsilon/sigmagsq;

div = 0.5*( sigmagsq/TotalNoise  -1  - log(sigmagsq/TotalNoise)   );


epsilonprime1 = 0;
epsilonprime2 = 2*log(2/hbinvdeltag);

for i = 1:max_iter

    epsilonprime = (epsilonprime1 + epsilonprime2)/2;
    
    ntemp =   log(2/hbinvdeltag) /  (   0.5*epsilonprime  - 0.5*log(1+epsilonprime)  );

    if logpe > log(hbinvdeltag/2) - 0.5*ntemp*epsilonprime*(sigmagsq/TotalNoise-1) - ntemp*div % The pe condition is satisfied; n can be made larger
    
        epsilonprime1 = epsilonprime;

    else
    
        epsilonprime2 = epsilonprime;

    end
    
end

n = - max(ntemp,1);

%REmember to return negative value for n