% MoM reference for the gprMax triangular bow-tie antenna.

clear;
close all;

script_dir = fileparts(mfilename("fullpath"));
results_dir = fullfile(script_dir, "results");
if ~isfolder(results_dir)
    mkdir(results_dir);
end

pattern_frequency = 0.82e9;
reference_impedance = 50;
outer_length = 101e-3;
outer_width = 100e-3;
flare_angle = 2 * atand(outer_width / outer_length);

% MATLAB's default bow tie lies in the y-z plane with its long axis along z.
% Rotate it about y to reproduce the x-y gprMax triangle pair.
antenna = bowtieTriangular( ...
    Length=outer_length, ...
    FlareAngle=flare_angle, ...
    Tilt=90, ...
    TiltAxis=[0 1 0] ...
);

set(groot, "defaultFigureVisible", "off");
mesh(antenna, MaxEdgeLength=4e-3);

alpha = (-180:2:180).';
theta = abs(alpha);
elevation = (-90:2:90).';
azimuth_columns = [0 180];

xz_grid = pattern( ...
    antenna, ...
    pattern_frequency, ...
    azimuth_columns, ...
    elevation, ...
    Type="directivity", ...
    Normalize=false ...
);
if size(xz_grid, 1) ~= numel(elevation)
    xz_grid = xz_grid.';
end
if ~isequal(size(xz_grid), [numel(elevation), numel(azimuth_columns)])
    error("Unexpected Antenna Toolbox elevation-pattern dimensions.");
end

cut_elevation = 90 - theta;
elevation_index = round((cut_elevation + 90) / 2) + 1;
xz_column = ones(size(alpha));
xz_column(alpha < 0) = 2;
xz_index = sub2ind(size(xz_grid), elevation_index, xz_column);
xz_directivity_dbi = xz_grid(xz_index);

xy_directivity_dbi = pattern( ...
    antenna, ...
    pattern_frequency, ...
    alpha, ...
    0, ...
    Type="directivity", ...
    Normalize=false ...
);
xy_directivity_dbi = reshape(xy_directivity_dbi, [], 1);
if numel(xy_directivity_dbi) ~= numel(alpha)
    error("Unexpected Antenna Toolbox azimuth-pattern dimensions.");
end

pattern_output = table( ...
    alpha, ...
    theta, ...
    cut_elevation, ...
    xz_directivity_dbi, ...
    xy_directivity_dbi, ...
    VariableNames=[ ...
        "angle_deg", ...
        "theta_deg", ...
        "elevation_deg", ...
        "xz_directivity_dbi", ...
        "xy_directivity_dbi" ...
    ] ...
);
writetable( ...
    pattern_output, ...
    fullfile(results_dir, "bowtie_antenna_matlab_pattern.csv") ...
);

gprmax_s11_path = fullfile(results_dir, "bowtie_antenna_gprmax_s11.csv");
if ~isfile(gprmax_s11_path)
    error( ...
        "Run bowtie_antenna_gprmax.py first so the shared frequency grid " + ...
        "is available at %s", ...
        gprmax_s11_path ...
    );
end
gprmax_s11 = readtable(gprmax_s11_path);
frequency_stride = 8;
frequency_indices = (1:frequency_stride:height(gprmax_s11)).';
if frequency_indices(end) ~= height(gprmax_s11)
    frequency_indices(end + 1) = height(gprmax_s11);
end
s11_frequency = gprmax_s11.frequency_hz(frequency_indices);
input_impedance = reshape(impedance(antenna, s11_frequency), [], 1);
s11 = (input_impedance - reference_impedance) ./ ...
    (input_impedance + reference_impedance);
s11_magnitude_db = 20 * log10(abs(s11));
s11_phase_deg = rad2deg(angle(s11));

s11_output = table( ...
    s11_frequency, ...
    real(input_impedance), ...
    imag(input_impedance), ...
    real(s11), ...
    imag(s11), ...
    s11_magnitude_db, ...
    s11_phase_deg, ...
    -s11_magnitude_db, ...
    VariableNames=[ ...
        "frequency_hz", ...
        "input_impedance_real_ohm", ...
        "input_impedance_imag_ohm", ...
        "s11_real", ...
        "s11_imag", ...
        "s11_magnitude_db", ...
        "s11_phase_deg", ...
        "return_loss_db" ...
    ] ...
);
writetable(s11_output, fullfile(results_dir, "bowtie_antenna_matlab_s11.csv"));

save( ...
    fullfile(results_dir, "bowtie_antenna_matlab.mat"), ...
    "pattern_frequency", ...
    "reference_impedance", ...
    "outer_length", ...
    "outer_width", ...
    "flare_angle", ...
    "alpha", ...
    "xz_directivity_dbi", ...
    "xy_directivity_dbi", ...
    "s11_frequency", ...
    "frequency_stride", ...
    "input_impedance", ...
    "s11" ...
);

fprintf( ...
    "Saved MATLAB bow-tie patterns at %.3f GHz and S11 at %d frequencies to %s\n", ...
    pattern_frequency / 1e9, ...
    numel(s11_frequency), ...
    results_dir ...
);