%% Figure 5.10 reproduction
%
clear;
beta = 50;
Nu = [1:10,20:10:100,200:100:1000];  % Number of Users
M = 32;     %Number of PPM positions
Ep = 1;     % Energy of pulse
No = 0.05;  % Noise value
simnum = 1000;
%
% rho
%
rho0=2*Ep/No;
%rho0=12.5;
rhoI=(3*beta)./((Nu-1)+(3*beta/rho0))  % channel symbol SNR
%
% r.v. simulation
%
total = 0;
for s=1:simnum
    sum =1;
    v1=normrnd(sqrt(rhoI),1);
    %v1 = 3.46;
    for m=2:M   
        vm=normrnd(0,1);
        %vm = 0;
        temp = exp(sqrt(rhoI).*(vm-v1));
        sum = sum + temp;
    end;
    CM = log2(M) - log2(sum);
    total = total + CM;
end
capacity = total/simnum
%
% plot the graph
%
semilogx(Nu,capacity,'-^')
grid;
ylabel('Capacity (bits/symbol)')
xlabel('Number of Users')
title('Number of Users vs Capacity (bits/symbol)')
hold;



M = 16;     %Number of PPM positions
%
% r.v. simulation
%
total = 0;
for s=1:simnum
    sum =1;
    v1=normrnd(sqrt(rhoI),1);
    for m=2:M   
        vm=normrnd(0,1);
        temp = exp(sqrt(rhoI).*(vm-v1));
        sum = sum + temp;
    end;
    CM = log2(M) - log2(sum);
    total = total + CM;
end
capacity = total/simnum
semilogx(Nu,capacity,'-v') % plot the graph for M=16




M = 8;     %Number of PPM positions
%
% r.v. simulation
%
total = 0;
for s=1:simnum
    sum =1;
    v1=normrnd(sqrt(rhoI),1);
    for m=2:M   
        vm=normrnd(0,1);
        temp = exp(sqrt(rhoI).*(vm-v1));
        sum = sum + temp;
    end;
    CM = log2(M) - log2(sum);
    total = total + CM;
end
capacity = total/simnum
semilogx(Nu,capacity,'-s') % plot the graph for M=8




M = 4;     %Number of PPM positions
%
% r.v. simulation
%
total = 0;
for s=1:simnum
    sum =1;
    v1=normrnd(sqrt(rhoI),1);
    for m=2:M   
        vm=normrnd(0,1);
        temp = exp(sqrt(rhoI).*(vm-v1));
        sum = sum + temp;
    end;
    CM = log2(M) - log2(sum);
    total = total + CM;
end
capacity = total/simnum
semilogx(Nu,capacity,'-+') % plot the graph for M=4






M = 2;     %Number of PPM positions
%
% r.v. simulation
%
total = 0;
for s=1:simnum
    sum =1;
    v1=normrnd(sqrt(rhoI),1);
    for m=2:M   
        vm=normrnd(0,1);
        temp = exp(sqrt(rhoI).*(vm-v1));
        sum = sum + temp;
    end;
    CM = log2(M) - log2(sum);
    total = total + CM;
end
capacity = total/simnum
semilogx(Nu,capacity,'-*') % plot the graph for M=2
hold off;