% Simalysis for bipolar network
% M. Haenggi, Oct. 2020
fs=18;
cols=[1 0 0;0 0 1;.49 .18 .56; .2 .7 .2];
N=50; % mean nr of points simulated
RUNS=500;
ALPHA=[2.1 2.25 2.5 3];
DELTA=2./ALPHA;
THETA=logspace(-2,2,50);
csp=zeros(length(ALPHA),length(THETA),RUNS);
cr2=N/pi; % squared simulation radius
meani=pi*DELTA./(1-DELTA).*cr2.^(1-1./DELTA); % mean interference not simulated
r=1/4; % link distance
a=0;
for k=1:RUNS
    n=poissrnd(N);
    x2=N/pi*rand(1,n); % squared distances (density 1)
    a=0;
    for alpha=ALPHA
      a=a+1;
      ra=r^alpha;
      i=0;
      for theta=THETA
          i=i+1;
          % double simalysis: avgeraging over Rayleigh fading and adding
          % the correction factor to account for interference outside
          % simulation region
          csp(a,i,k)=1/prod(1+theta*ra*x2.^(-1/DELTA(a)))*exp(-theta*ra*meani(a));
      end
   end
end
sp=mean(csp,3);
T=[0 THETA./(1+THETA) 1]; % MH scale
a=0;
leg=[];
h=0*ALPHA; g=h;
for alpha=ALPHA
    a=a+1;
    delta=2/alpha;
    figure(1)
    h(a)=plot(10*log10(THETA),sp(a,:),'LineWidth',2,'Color',cols(a,:));
    hold on
    ana=exp(-pi*THETA.^delta*r^2*gamma(1+delta)*gamma(1-delta));
    plot(10*log10(THETA),ana,'--','Color',cols(a,:),'LineWidth',2);
    figure(2)
    g(a)=plot(T,[1 sp(a,:) 0],'LineWidth',2,'Color',cols(a,:));
    hold on
    plot(T,[1 ana 0],'--','Color',cols(a,:),'LineWidth',2);
    leg=[leg '\alpha=' num2str(alpha) ' '];
end
leg=leg(1:end-1);
set(gca,'FontSize',fs);
xlabel('\theta [MH]','FontSize',fs);
ylabel('$\bar F_{\rm{SIR}}(\theta)$','FontSize',fs,'interpreter','latex');
grid on
hold off
legend(g,split(leg),'FontSize',fs);
figure(1)
set(gca,'FontSize',fs);
xlabel('\theta [dB]','FontSize',fs);
ylabel('$\bar F_{\rm{SIR}}(\theta)$','FontSize',fs,'interpreter','latex');
grid on
hold off
legend(h,split(leg),'FontSize',fs);