clear
Nsample = 30000;
N = 3; 
g = 6;
alpha = 4/g;
lamTotal = 0.001;
Threshold = 2*ones(1,N);
P = 1;M = 64;
A = 1.0;
sigG = 1.0;
sigH = 1.0;
xVec = linspace(0,lamTotal,11);

for m = 2:8
    N=m
    lamC = lamTotal/20;
    lamI = lamTotal - lamC;
    
    C0 = 4*(A*sigG/(sqrt(P)*sigH))^(-alpha)*(1/(cos(pi*alpha/4)*sqrt(pi)))*gamma(1/2+alpha/2)*gamma(1-alpha/2);
    
    %SC Receiver
    temp1 = 0;
    for count = 1:N
        temp1 = temp1 + (-1)^(count+1) * nchoosek(N,count)*count^(alpha/2);
    end
    C_SC = temp1*C0;
    outageSC_Theo(m) = C_SC*(lamI + lamC)*Threshold(1)^(alpha/2);
    
    %EGC Receiver
    outageEGC_Theo(m) = C0*(N^(1-alpha/2)*lamI + lamC)*Threshold(1)^(alpha/2);
       
    %MRC Receiver
    h = randn([Nsample N]) + 1i*randn([Nsample N]);
    Iconst = mean( sum( abs(h).^(alpha),2 )./ (sum( abs(h).^(2),2 ).^(alpha) ) );
    Cconst = mean( 1./sum( abs(h).^(2),2 ).^(alpha/2) );
    Cmrc = 4*(A*sigG/(sqrt(P)))^(-alpha)*(1/(cos(pi*alpha/4)))*gamma(1/2+alpha/2)*2^(alpha/2)/sqrt(pi);
    outageMRC_Theo(m) = Cmrc*(lamI*Iconst + lamC*Cconst)*Threshold(1)^(alpha/2);
 
    %PDC Receiver
    tempC = 0;
    for k = 1:N
        tempC = tempC + (-1)^(k+1)*nchoosek(N,k)*prod( (1:k-1) + alpha/2)/factorial(k-1);
    end
    outagePDC_Theo(m) = 0;
    outagePDC_Theo(m) =  outagePDC_Theo(m) + tempC*(A*sigG/(sqrt(P)*sigH))^(-alpha)*alpha*pi^2/(2*sin(alpha*pi/2))*lamC*Threshold(1)^(alpha/2);
    outagePDC_Theo(m) =  outagePDC_Theo(m)+ (A*sigG/(sqrt(P)*sigH))^(-alpha)*(alpha*(pi^2)/(2*sin(alpha*pi/2)))^N *(lamI)^N *Threshold(1)^(N*alpha/2);

end

load OutageVsNmix_Sim;
outageEGC = outageSim(:,2)';
outageSC = outageSim(:,3)';
outageMRC = outageSim(:,4)';
outagePDC = outageSim(:,5)';
outageHL = outageSim(:,6)';
outageSL = outageSim(:,7)';

close all
figure
hold on
plot(log10(outageSC),'bo');
plot(log10(outageSC_Theo));
plot(log10(outageMRC),'r*');
plot(log10(outageMRC_Theo),'r-');
plot(log10(outageEGC),'gs');
plot(log10(outageEGC_Theo),'g-');
plot(log10(outagePDC),'kd');
plot(log10(outagePDC_Theo),'k-');
plot(log10(outageHL),'m');
plot(log10(outageSL),'c');

outageAllNmix = [(1:8)' outageEGC' outageEGC_Theo' outageSC' outageSC_Theo' outageMRC' outageMRC_Theo' outagePDC' outagePDC_Theo' outageHL' outageSL']

save OutageVsNmixAll outageAllNmix