close all ;
clear; clc;

%% Parametry systemu
frequencies = [865e6, 915e6];               % [Hz], podstawowe cz?stotliwo?ci
frequency_range = 800e6:5e6:1000e6;         % [Hz], zakres modelowania
rPwrMax = 5;                                % Maksymalny promie? [m]
r_values = [0.9, 0.7, 0.5] * rPwrMax;       % Promienie analizy
mode_tag_methods = {'real', 'ideal'};      % Dwie metody modelowania anteny taga

%% Parametry anteny czytnika
G0R_dBi = 7;                                % Maksymalny zysk [dBi]
HPBWR = 58;                                 % Half-Power Beamwidth [deg]

%% Parametry ogólne
C_chip = 0.63e-12;  % [F] pojemno?? chipa RFID
Rc = 150;        % [Ohm] rezystancja chipa
Za = 50;         % [Ohm] impedancja anteny (przyj?ta sta?a)
RC = 0.5;                                   % strata polaryzacyjna [dB]
Ptx_dBm = 30;                               % Moc nadajnika [dBm]
Pr_min_dBm = -21;                           % Czułość taga [dBm]

%% Wczytanie danych z 4nec2 (charakterystyka promieniowania taga)

rawData = readtable('book3.csv', 'Delimiter', ';');
theta = rawData{:,1};
phi = rawData{:,2};
gain = rawData{:,3};
gain(gain < -100) = NaN;
valid = ~isnan(gain);
theta = theta(valid);
phi = phi(valid);
gain = gain(valid);
[uniquePairs, ia, ~] = unique([theta, phi], 'rows');
theta = uniquePairs(:,1);
phi = uniquePairs(:,2);
gain = gain(ia);
F_tag = scatteredInterpolant(theta, phi, gain, 'linear', 'nearest');

%% Funkcja anteny czytnika (modelowana)
readerPattern = @(theta, G0, HPBW) G0 - 3*(abs(theta)/HPBW).^(-log(0.5)/log(cosd(HPBWR/2)));
idealTagPattern = @(theta) 0;  % Model izotropowy (sta?y zysk 0 dBi)

%% Analiza i wykresy: PT(f) dla ró?nych r i metod

for m = 1:length(mode_tag_methods)
    method = mode_tag_methods{m};
    figure('Name', sprintf('Metoda: %s', method)); hold on;
    legends = {};

    for r = r_values
        PT = zeros(size(frequency_range));
        for i = 1:length(frequency_range)
            f = frequency_range(i);
            lambda = 3e8 / f;
            gain_reader = readerPattern(30, G0R_dBi, HPBWR);

            if strcmp(method, 'real')
                gain_tag = F_tag(0, -30);
            else
                gain_tag = idealTagPattern(0);
            end

            L = (4 * pi * r / lambda)^2;
            L_dB = 10 * log10(L);
          %% Obliczenie dopasowania impedancyjnego
Xc = 2 * pi * f * C_chip;
Zc = Rc + 1i * Xc;
Gamma = (Zc - Za) / (Zc + Za);
matching_loss = 1 - abs(Gamma)^2;

% Moc dostarczona z uwzgl?dnieniem odbicia
Pr_lin = 10^((Ptx_dBm + gain_reader + gain_tag - L_dB - RC)/10);
Pd_lin = Pr_lin * matching_loss;
PT(i) = 10 * log10(Pd_lin);
        end
        plot(frequency_range/1e6, PT);
        legends{end+1} = sprintf('r = %.1f m', r);
    end

    xlabel('Cz?stotliwo?? [MHz]');
    ylabel('Moc dostarczona do taga [dBm]');
    title(sprintf('Charakterystyka PT(f) – metoda: %s', method));
    legend(legends, 'Location', 'best');
    grid on;
end

%% Mapa 2D mocy odbieranej Pr(x, y) dla f = 865 MHz
f_map = 865e6;
lambda = 3e8 / f_map;
range = -10:0.2:15;
[X, Y] = meshgrid(range, range);

for m = 1:length(mode_tag_methods)
    method = mode_tag_methods{m};
    Pr_map = zeros(size(X));

    for i = 1:size(X, 1)
        for j = 1:size(X, 2)
            x = X(i, j);
            y = Y(i, j);
            r = sqrt(x^2 + y^2);
            theta = atan2d(y, sqrt(x^2));
            theta = mod(theta + 360, 360);

            gain_reader = readerPattern(theta, G0R_dBi, HPBWR);
            if strcmp(method, 'real')
                gain_tag = F_tag(theta, -30);
            else
                gain_tag = idealTagPattern(theta);
            end

            L = (4 * pi * r / lambda)^2;
            L_dB = 10 * log10(L);
          %% Impedancyjne dopasowanie
Xc = 2 * pi * f_map * C_chip;
Zc = Rc + 1i * Xc;
Gamma = (Zc - Za) / (Zc + Za);
matching_loss = 1 - abs(Gamma)^2;

% Moc z uwzgl?dnieniem odbicia
Pr_lin = 10^((Ptx_dBm + gain_reader + gain_tag - L_dB - RC)/10);
Pd_lin = Pr_lin * matching_loss;
Pr_map(i, j) = 10 * log10(Pd_lin);
        end
    end

    figure;
    imagesc(range, range, Pr_map .* (Y >= 0));
    colormap('jet');         % lub 'parula', 'hot', 'turbo'
    caxis([-40 0]);          % ustalony zakres dBm
    hold on;
    [C, h] = contour(X, Y, Pr_map, [-21 -21], 'LineColor', 'k', 'LineWidth', 1.5);
    clabel(C, h, '-18 dBm');
    set(gca, 'YDir', 'normal');
    colorbar;
    title(sprintf('Mapa Pr(x,y) dla f = %.0f MHz, metoda: %s', f_map/1e6, method));
    xlabel('x [m]'); ylabel('y [m]');
end

%% Ocena poprawnej pracy dla f = 865, 915 MHz
theta_tag = 90;
phi_tag = -30;

fprintf('\n--- Ocena poprawnej pracy (f = 865, 915 MHz) ---\n');
for f = frequencies
    lambda = 3e8 / f;
    fprintf('\nCzestotliwosc: %.0f MHz\n', f/1e6);
    for r = r_values
        gain_reader = readerPattern(theta_tag, G0R_dBi, HPBWR);
        gain_tag = F_tag(theta_tag, phi_tag);
        L = (4 * pi * r / lambda)^2;
        L_dB = 10 * log10(L);
       %% Dopasowanie impedancji
Xc = 2 * pi * f * C_chip;
Zc = Rc + 1i * Xc;
Gamma = (Zc - Za) / (Zc + Za);
matching_loss = 1 - abs(Gamma)^2;
Pr_lin = 10^((Ptx_dBm + gain_reader + gain_tag - L_dB - RC)/10);
Pd_lin = Pr_lin * matching_loss;
Pr_dBm = 10 * log10(Pd_lin);
        status = 'NIE';
        if Pr_dBm >= Pr_min_dBm
            status = 'TAK';
        end
        fprintf('r = %.2f m, Pr = %.2f dBm, Poprawna praca: ''%s''\n', r, Pr_dBm, status);
    end
end