%仿真代碼的初始化
clear;
r_cell=1000;
n_co_ch_users=6;
P_BS=0;
P_MS=0;
corr_fwd=0.0;
corr_rev=0.0;
K=1;
in_beam=0;
sector_min=zeros(6,3);
sector_max=zeros(6,3);
sector_min(:,1)=(pi/3)*[-3:2]';
sector_min([1:3],2)=pi/3*[-3 -1 1]';
sector_min(1,3)=-pi;
sector_max(:,1)=sector_min(:,1)+pi/3;
sector_max([1:3],2)=sector_min([1:3],2)+2*pi/3;
sector_max(1,3)=pi;
phi_center=zeros(6,3);
phi_center(:,1)=(pi/3)*[-3:2]'; 
phi_center([1:3],2)=(pi/3)*[-3 -1 1]';
phi_BW=[1 2 6]*pi/3;
num_sectors=[6 3 1];
num_snapshots=input('Number of snapshots= ');
cluster_size=input('Cluster size(3,4or7)=');
n_path=input('Path loss exponent=');
sigma_int=input('Shadowing std deviation-interference(dB)=');
sigma_S=input('Shadowing std deviation-desired signal(dB)=');
sec=input('Sectorization…(1=>60 degree,2=>120 degree,3=>omni).enter:');
ftb=input('Front-to-back ratio of the BS antennas (dB)=');
out_beam=in_beam-ftb;

%用直角坐標和極坐標表示仿真中基站和移動臺的位置
theta_N=[pi/6 0 pi/6 asin(1/(2*sqrt(7)))];
theta=pi/3*[0:5]';
aux_1=[1 0 2 3 0 0 4];
ind=aux_1(cluster_size);
bs_position=[sqrt(3*cluster_size)*r_cell*cos(theta+...
    theta_N(ind)) sqrt(3*cluster_size)*r_cell*sin(theta+...
    theta_N(ind))];

%把移動臺放入每一個小區中。假設移動臺在小區內均勻分布
sector=unidrnd(num_sectors(sec));
des_user_beta=rand(1)*phi_BW(sec)+phi_center(sector,sec);
des_user_r=sqrt(rand(1).*(r_cell^2));
co_ch_user_beta=rand(6,1)*phi_BW(sec)+phi_center(sector,sec);
co_ch_user_r=sqrt(rand(6,1)*(r_cell));

%計算同頻移動臺和所研究的移動臺的坐標位置
des_user_position=des_user_r*[cos(des_user_beta)...
    sin(des_user_beta)];
co_ch_user_position=[co_ch_user_r.*cos(co_ch_user_beta)...
    co_ch_user_r.*sin(co_ch_user_beta)]+bs_position;

%計算前向鏈路信號的各階距
m_S_fwd=P_BS-10*K*n_path*log10(des_user_r);
aux_01=((des_user_position(1)-bs_position(:,1))+...
    sqrt(-1)*(des_user_position(2)-bs_position(:,2)));
beta_fwd=angle(aux_01);
d_I_fwd=abs(aux_01);
clear gain_fwd
for k=1:n_co_ch_users
if(beta_fwd(k)>=...
        sector_min(sector,sec))&(beta_fwd(k)<...
        sector_max(sector,sec))
gain_fwd(k)=in_beam;
else
gain_fwd(k)=out_beam
end
end
m_I_fwd=P_BS-10*K*n_path*log10(d_I_fwd)+gain_fwd.';
sigma_I_fwd=sigma_int*ones(length(m_I_fwd),1);
[m_I_total_fwd,sigma_I_total_fwd]=...
    wilkinson(m_I_fwd,sigma_I_fwd,corr_fwd);

%計算反向鏈路信號的各階距
m_S_rev=P_MS-10*K*n_path*log10(des_user_r);
aux_02=(co_ch_user_position(:,1)+...
    sqrt(-1)*co_ch_user_position(:,2));
beta_rev=angle(aux_02);
d_I_rev=abs(aux_02);
clear gain_rev
for k=1:n_co_ch_users
if(beta_rev(k)>=...
sector_min(sector,sec))&(beta_rev(k)<...
sector_max(sector,sec))
gain_rev(k)=in_beam;
else
gain_rev(k)=out_beam;
end
end
m_I_rev=P_MS-10*K*n_path*log10(d_I_rev)+gain_rev.';
sigma_I_rev=sigma_int*ones(length(m_I_rev),1);
[m_I_total_rev,sigma_I_total_rev]=...
    wilkinson(m_I_rev,sigma_I_rev,corr_rev);

%計算SIR統計特性
for i=1:1000

m_SIR_fwd(i)=m_S_fwd-m_I_total_fwd;
sigma_SIR_fwd(i)=sqrt(sigma_S^2+sigma_I_total_fwd^2-...
    2*corr_fwd*sigma_S*sigma_I_total_fwd);
m_SIR_rev(i)=m_S_rev-m_I_total_rev;
sigma_SIR_rev(i)=sqrt(sigma_S^2+sigma_I_total_rev^2-...
    2*corr_rev*sigma_S*sigma_I_total_rev);

%中斷概率
inside='(1/(sigma_SIR_fwd(i)*sqrt(2*pi))*exp(-(x-m_SIR_fwd(i))^2)/(2*sigma_SIR_fwd(i)^2))';
P_outage=int(inside,0,10); 

end

% 畫圖
C_1=[1:10];
A_1=linspace(0.1,10,50);
C_2=[12:2:20];
A_2=linspace(3,20,50);
C_3=[30:10:100];
A_3=linspace(13,120,50);
for i=1:length(C_1)
    for j=1:length(A_1)
        p_1(j,i)=erlang_b(A_1(j),C_1(i));
    end
end
for i=1:length(C_2)
    for j=1:length(A_2)
        p_2(j,i)=erlang_b(A_2(j),C_2(i));
    end
end
for i=1:length(C_3)
    for j=1:length(A_3)
        p_3(j,i)=erlang_b(A_3(j),C_3(i));
    end
end
x1=[.1 .1 .2 .2 .3 .3 .4 .4 .5 .5 .6 .6 .7 .7 .8 .8 .9 .9];
y1=10.^[-4 1 1 -4 -4 1 1 -4 -4 1 1 -4 -4 1 1 -4 -4 -1];
y2=[.1 .1 .2 .2 .3 .3 .4 .4 .5 .5 .6 .6 .7 .7 .8 .8 .9 .9];
x2=10.^[-1 3 3 -1 -1 3 3 -1 -1 3 3 -1 -1 3 3 -1 -1 3];
loglog(A_1,p_1,'k-',A_2,p_2,'k-',A_3,p_3,'k-',...
    x1,y1,'k:',10*x1,y1,'k:',...
    100*x1,y1,'k:',1000*x1,y1,'k:',...
    x2,y2,'k:',x2,0.1*y2,'k:',x2,0.01*y2,'k:');
xlabel('Offered traffic (Erlangs)')
ylabel('Blocking probability')
axis([0.1 120 0.001 0.1])
text(.115, .115,'C=1')
text(.6, .115,'C=2')
text(1.18, .115,'3')
text(2, .115,'4')
text(2.8, .115,'5')
text(7, .115,'10')
text(9, .115,'12')
text(17, .115,'20')
text(27, .115,'30')
text(45, .115,'50')
text(100, .115,'100')