wcn_code.m

clc;
clear;
close all;
radius = 4700000/2; %m
center_x(1) = 0;
center_y(1) = 0;
%FORMING THE CELLS
n=7;               %number of cells
D=radius*sqrt(3);  %distance btw cells
t = (0:pi/3:2*pi);
%center cell
x1 = center_x(1) + radius*cos(t);
y1 = center_y(1) + radius*sin(t);
plot(x1,y1)
hold on
%adjacent cells
N=5;                             %the second cell corresponds to the upper one
directions = (30:60:360);
adj_center_x=zeros([1 5]);
adj_center_y=zeros([1 5]);
for i=2:(N+1)
    center_x(i) = center_x(1)+D*cosd(directions(i));
    center_y(i) = center_y(1)+D*sind(directions(i));
    x = center_x(i) + radius*cos(t);
    y = center_y(i) + radius*sin(t);
    plot(x,y)
    adj_center_x(1,i-1)=center_x(i);
    adj_center_y(1,i-1)=center_y(i);
end
hold on;
%Right Upper cell
x2= (radius+(radius/2)) + radius*cos(t);
y2 = ((sqrt(3)/2)*radius) + (radius*sin(t));
cx4 = (radius+(radius/2));  %BS Location Right Upper cell x-axis
cy4 = ((sqrt(3)/2)*radius);  %BS Location Right Upper cell y-axis
plot(x2,y2);
hold on;
%upper additional cells
N=5;
directions = (30:30:360);
adj_center_x2=zeros([1 5]);
adj_center_y2=zeros([1 5]);
for i=2:(N+1)
    if i==2 || i==4 || i==6
    center_x(i) = center_x(1)+sqrt(3)*D*cosd(directions(i));
    center_y(i) = center_y(1)+sqrt(3)*D*sind(directions(i));
    x = center_x(i) + radius*cos(t);
    y = center_y(i) + radius*sin(t);
    else
    center_x(i) = center_x(1)+2*D*cosd(directions(i));
    center_y(i) = center_y(1)+2*D*sind(directions(i));
    x = center_x(i) + radius*cos(t);
    y = center_y(i) + radius*sin(t);
    end
adj_center_x2(1,i-1)=center_x(i);
adj_center_y2(1,i-1)=center_y(i);
plot(x,y)
end


%BS DETERMINATION
cx1=0;
cy1=0;  %center cell
%second BS is the upper one then it goes counter clockwise
%the right upper cell (cx4,cy4) was added later
%the 8th BS corresponds to the right most upper one then it goes counter
%clockwise
bs_x = [cx1 adj_center_x cx4 adj_center_x2];
bs_y = [cy1 adj_center_y cy4 adj_center_y2];
hold on;
plot(bs_x,bs_y,'r+')
BS = [bs_x' bs_y'];
%numbering the cells
sayi = [1,2,3,4,5,6,7,8,9,10,11,12];
for i=1:12     % enumeration of the cells
    text(bs_x(i)+75100,bs_y(i)-75100,(num2str(sayi(i))))
end


%DISTRIBUTING USERS
user = user_distribution(1,BS,radius-1000000);
user_x = user(:,1);
user_y = user(:,2);
hold on;
plot(user_x,user_y,'b*')
user_x_initial = user_x;
user_y_initial = user_y;


%ADDING MOBILITY & DISTANCE CALCULATION
distan = [12;12;35];
coeff = 3*(rand(length(user_x),1) - 0.5);
coeff2 = 4*(rand(length(user_x),1) - 0.5);
start = user_y - coeff.*user_x;
R=1e7;
m = 1;
plot(user_x,user_y,'o')  %ploting points old points
axis([-1.5*R 1.5*R -1.5*R 1.5*R])  %limiting windomw
drawnow
while m<=35
    user_x =user_x+ R/100.*step_func(coeff);  %.*(-1).^ceil(2*rand(1,n));% adding random velocity. Factor of 20 was picked arbitrarily
    user_y = coeff.*user_x + start;  %.*(-1).^ceil(2*rand(1,n));% adding random velocity. BURAYI EKLEME DENE coeff2.*user_x.^2 +
    plot(user_x,user_y,'x')   %ploting new points
    axis([-1.5*R 1.5*R -1.5*R 1.5*R])
    drawnow  %display change now    
    user = [user_x user_y];
    for i=1:length(user_x)
        for j=1:length(bs_x)
            a = sqrt((user_x(i) - bs_x(j))^2 + (user_y(i) - bs_y(j))^2);
            distan(j,i,m) = a;    
        end
    end   
    m = m+1;
end


%PATH LOSS CALCULATION
h=781000;  %satellite height in m
h2= 8534;  %UAV height in m
delta_h = h-h2;
freq = 1620;  %Mhz
%free space loss in dB
d_up=zeros(12,12);
free_space=zeros(12,12);
L=zeros(12,12);
for i=1:12
    for j = 1:12
        for k=1:35
            d_up(i,j,k) = sqrt(distan(i,j,k)^2+delta_h^2)*10^-3;
            free_space(i,j,k) = 32.45+20*log10(freq)+20*log10(d_up(i,j,k)); 
            %basic path loss = FSPL + SF  (SF is random number generated by the normal distribution)
            %elevation angle of UAV is 12 deg    
            pd = makedist('Normal','sigma',1.79);   %for suburban,elevation=10 sigma=1.79 (dB)
            %rng('default')   
            sf=random(pd);
            L(i,j,k) = free_space(i,j,k) + sf;  
            
            %received power Pr = Pt+Gt+Gr-FSPL-Other losses
            losses_up = 0.8;  %dB
            uav_gain = 3; %dBi
            satellite_antenna_rx = 25.2; %dBi
            uav_power = 20; %dBW
            Pr_up(i,j,k) = uav_power+uav_gain+satellite_antenna_rx-free_space(i,j,k)-losses_up;
             
            losses_down=0.8; %dB
            satellite_antenna_tx = 41; %dBi
            satellite_tx_power = -29.18; %dBW
            Pr_down(i,j,k) =satellite_tx_power+uav_gain+satellite_antenna_tx-free_space(i,j,k)-losses_down;
            
        end
    end    
end


%SNR CALCULATION
NF = 7 ; %noise figure in dB
bw = 5;  %MHz
n0_downlink = -205.93;  %dBm/Hz 
n0_uplink = -202.3; %dBm/Hz
Pn_down = NF + n0_downlink + 10*log10(bw*10^6); %receiver noise power at UE in dBm
Pn_up = NF + n0_uplink + 10*log10(bw*10^6);
SNR_downlink = Pr_down - Pn_down;  %snr of received signal in dB
SNR_uplink = Pr_up- Pn_up;



%RSS CALCULATION
for i=1:12
    for j=1:12
        for k=1:35
            received_signal = -104.27;  %dBW
            rss(i,j,k) = received_signal - 20*log10(d_up(i,j,k))+satellite_antenna_tx+3;
        end
    end
end
%each user rss according to 12 BS in 30 iterations
user1=squeeze(rss(:,1,:));user2=squeeze(rss(:,2,:));
user3=squeeze(rss(:,3,:));user4=squeeze(rss(:,4,:));
user5=squeeze(rss(:,5,:));user6=squeeze(rss(:,6,:));
user7=squeeze(rss(:,7,:));user8=squeeze(rss(:,8,:));
user9=squeeze(rss(:,9,:));user10=squeeze(rss(:,10,:));
user11=squeeze(rss(:,11,:));user12=squeeze(rss(:,12,:));
user_rss = [user1;user2;user3;user4;user5;user6;user7;user8;user9;user10;user11;user12];
user_rss = mat2cell(user_rss,[12 12 12 12 12 12 12 12 12 12 12 12]);
%plotting the rss values
for m=1:12
   figure(m+1);
   plot(1:35,squeeze(rss(:,m,:)))
   ylabel('RSS')
   title(sprintf('%dth user RSS values',m))
   legend({'BS1','BS2','BS3','BS4','BS5','BS6','BS7','BS8','BS9','BS10','BS11','BS12'},'Location','northwest','NumColumns',6)
end



%HANDOVER DECISION FOR EACH USER USING RSS VALUES AND DISTANCE
position_min_distance = zeros(35,12);
min_distance = zeros(35,12);
for k=1:35
    for l=1:12
        min_distance(k,l) = min(distan(:,l,k));
        position_min_distance(k,l) = find(distan(:,l,k) == min_distance(k,l));
    end 
end

current_base_station = zeros(12,1);
maximum = zeros(12,1);
for i = 1:12
     maximum(i) = max(rss(:,i,1));
     if find(rss(:,i,1) == maximum(i)) ==  min_distance(i,1)
        current_base_station(i) =  min_distance(i,1);
     else
         current_base_station(i) = find(rss(:,i,1) == maximum(i));
     end
end

margin = 5;
handover_rate = 0;
bs_no=ones(12,35);
finishNow = false;
difference=zeros(35,12);
ttt=5;
ttt_initial = zeros(1,12);

for m=1:12
    finishNow = false;
    ho=zeros(1,5);
     for k=1:30 
             user_m = user_rss(m,:);          %users are selected
             user_m = cell2mat(user_m);       %in each row 12 BS are represented
             max_user_m = max(user_m(:,k));
             difference(k,m) = max_user_m - user_m(current_base_station(m),k);
             if difference(k,m) > margin && position_min_distance(k,m) == find(user_m(:,k) == max_user_m  )
                 ttt_initial(m) = k;  
                 for l=1:ttt
                     difference2(1,l)= max(user_m(:,ttt_initial(m)+l)) - user_m(current_base_station(m),ttt_initial(m)+l);
                     if difference2(1,:) > margin 
                         ho(l) = 1;
                     end
                     if prod(ho,2) == 1
                         handover_rate = handover_rate+1;
                         finishNow=true;
                     end
                 end
                  current_base_station(m) =  find(user_m(:,k) == max_user_m);
             end 
         if finishNow
             break
         end
     end
end


ttt_initial
handover_rate


%% used functions
function [user] = user_distribution(user_no_perbs,bs,cell_rad)
        
k=1;
for i=1:length(bs)
    for j=1:user_no_perbs
    [coord_x,coord_y] = get_coordinates(bs(i,1),bs(i,2), cell_rad);
    user(k,1)=coord_x;
    user(k,2)=coord_y;
    k=k+1;
    end
end
end

function [coord_x,coord_y] = get_coordinates(bs_x,bs_y,radius)
m1 = tan(30*pi/180);
m2 = tan(150*pi/180);
outside = 2;
R=radius;
while outside
     x_temp = R - 2*rand*R;
     y_temp = R - 2*rand*R;
     if (x_temp < sqrt(3)/2*R && x_temp>-sqrt(3)/2*R) && (y_temp < m1*x_temp+R && y_temp <m2*x_temp+R && y_temp > m2*x_temp-R && y_temp>m1*x_temp-R )
         outside = 0;
     end
end
coord_x = y_temp + bs_x;
coord_y = x_temp + bs_y;
end

function [ret] = step_func(array)
    
range = length(array);
  for i=1:range
      if array(i) > 0
          array(i) = 1;
      elseif array(i) < 0
          array(i) = -1;
      else 
          array(i) = 0;
      end
 
   end
ret = array;
end