clear all; clc; close all;

% figure 전체화면 setting
scrsz=get(0,'screensize');
fig1=figure('Position',[ 1 1 scrsz(3) scrsz(4)]);


theta = linspace(0, 2*pi);
radius = 1000; % 전체 지역 내 반경(m)
bound = 50; % tx와 rx가 존재하는 최대 반경(m)
tr_num = 100; % 시뮬레이션을 위한 tx+rx수
% turn_num = 10000; % h값을 랜덤하게 선택하기 위한 샘플링 횟수
delay_time = 0.0001;
beta = 2.5;
n0 = 10^-8; % channel의 noise 값
tx_point = zeros(1, 2, tr_num); % tx 위치 
rx_point = zeros(1, 2, tr_num); % rx 위치
d = zeros(tr_num, tr_num); % distance between devices
h = zeros(tr_num, tr_num); % channel
p_sample = 100; % sampling num
p = linspace(0.002, 0.5, p_sample); % power[W]

X = radius*cos(theta);
Y = radius*sin(theta);
subplot(2, 2, 2);
subplot(2, 2, 3);
subplot(2, 2, 4);
subplot(2, 2, 1);
plot(X, Y, 'b--');
title("Simulation region");
xlabel("x");
ylabel("y");
hold on;

h = sqrt(1/2) * randn(tr_num, tr_num) + i*randn(tr_num, tr_num);


% 랜덤 좌표로 tx와 rx위치를 찍어준다.
for j = 1:tr_num
    while 1
        rand_rad = (radius - bound) * rand;
        rand_theta = 2 * pi * rand;
    
        txx = rand_rad*cos(rand_theta);
        txy = rand_rad*sin(rand_theta);
        flag = 0;
        
        k = 1;
        while k < j
                tar = tx_point(:, :, k);
            
            tar_x = tar(1, 1);
            tar_y = tar(1, 2);

            det = (txx-tar_x)^2 + (txy-tar_y)^2 - (2*bound)^2;
            
            k = k + 1;
            
            if det <= 0 % bound가 겹치는 경우 다시 돌아야함
                flag = 1;
                break;
            else
                continue;
            end
        end
        
        if flag
            continue;
        else
            break; % flag 값이 0이면 바운드 내에 겹치는 점이 없다는 뜻이므로 break한뒤 좌표 저장
        end
    end
    
    tx_point(:, :, j) = [txx txy];
    
    pause(delay_time);
    
    subplot(2, 2, 1);
    plot(txx, txy, 'r*'); 
    
    % for indexing points
%     t = text(txx, txy, num2str(j));
%     t.Color = [1 0 0];
    %
    
    B1 = bound*cos(theta);
    B2 = bound*sin(theta);
    plot(B1+txx, B2+txy, 'b-');
    
    rand_rad = bound * rand;
    rand_theta = 2 * pi * rand;
    
    rxx = rand_rad*cos(rand_theta);
    rxy = rand_rad*sin(rand_theta);
    rx_point(:, :, j) = [rxx + txx, rxy + txy];
   
    subplot(2, 2, 1);
    plot(rxx + txx, rxy + txy, 'm^');
    
    
    
    for k = 1:j % tx
        txs = tx_point(1, 1, k);
        tys = tx_point(1, 2, k);    
        for w = 1:j % rx
            rxf = rx_point(1, 1, w);
            ryf = rx_point(1, 2, w);
            %calculate the distance between tx and rx
            d(k, w) = sqrt((txs-rxf)^2 + (tys-ryf)^2);
        end
    end
    
    htr = (d.^(-beta)).*abs(h).^2;
    mean_r = zeros(1, p_sample);
    
    for k = 1:j
        intf = n0;
        for rx = 1:j
            if k == rx
                continue;
            else
                intf = intf + htr(k, rx);
            end
        end
        r = log2(1+(htr(k, k)*p)/(intf+n0));
        mean_r = mean_r + r;
    end
    mean_r = mean_r/j;
    subplot(2, 2, 2);
    plot(p, mean_r, 'r-');
    title("전력 P에 따른 data rate 변화");
    xlabel("P");
    ylabel("rate");
    hold on;
    
    subplot(2, 2, 3);
    s1 = round(p_sample/3);
    plot(j, mean_r(s1), 'bo');
    title("단말장치 수에 따른 전력: " + p(s1)+"W 에서의 rate");
    xlabel("단말장치 수");
    ylabel("rate");
    hold on;
    
    plot([j j], [mean_r(s1), 0], 'b-');
    subplot(2, 2, 4);
    s2 = round(2*p_sample/3);
    plot(j, mean_r(s2), 'bo');
    title("단말장치 수에 따른 전력: " + p(s2)+"W 에서의 rate");
    xlabel("단말장치 수");
    ylabel("rate");
    hold on;
    plot([j j], [mean_r(s2), 0], 'b-');
end