% Vivaldi_Sweep.m
%
% Automated parameter sweep for the CPW-fed Vivaldi antenna, built on top
% of the now-debugged, fast-converging core from Vivaldi_OpenEMS_Start.m.
% Runs several (taper.L, taper.ya) combinations back-to-back with NO user
% interaction needed, and prints a summary table at the end so you can see
% at a glance which combination gets closest to S11 < -10 dB from 850 MHz
% upward - instead of testing one combination at a time by hand.
%
% Fixed/confirmed-stable choices carried over from Vivaldi_OpenEMS_Start.m:
%   - taper.y0 derived from cpw.center_w/cpw.gap (no short at the throat)
%   - taper.R = 0.023 (0.014/0.018 were numerically unstable)
%   - N = 10 taper polygon points, pre_tol/merge_tol mesh-line cleanup
%     (avoids the SmoothMesh hangs seen with denser point sets)
%   - coarser exploration mesh divisors (/20, /10)
% Each run's NrTs cap is deliberately smaller than the production script
% (1,000,000) to bound worst-case time per combination; if a combination
% is unstable it will hit this cap quickly rather than run for 30+ min.
%
% NOTHING in this script requires you to press anything or check anything
% mid-run. Start it, walk away, come back to the summary table at the end.
%
% Figures are suppressed during the sweep for speed; only the final S11
% curve for the automatically-selected best combination is plotted at the
% end (Figure 1).

close all
clear
clc

max_recursion_depth(10000);

physical_constants;
unit = 1e-3;

f_low  = 850e6;
f_high = 2000e6;
f0 = (f_low + f_high)/2;
fc = (f_high - f_low)/2 + 0.3e9;

substrate.epsR      = 4.3;
substrate.tand      = 0.023;
substrate.kappa     = substrate.tand * 2*pi*f0 * substrate.epsR * EPS0;
substrate.thickness = 1;

feed.throat_run = 15;
feed.R_ref      = 50;
cpw.center_w = 3.0;
cpw.gap      = 0.4;
taper.y0 = cpw.center_w/2 + cpw.gap;
taper.R  = 0.023; % confirmed stable - not swept

% ---- combinations to test: [taper.L, taper.ya] in mm ----
combos = [ ...
    160  50; ...
    200  50; ...
    240  50; ...
    160  60; ...
    200  60; ...
    240  60  ...
];

N = 10;
NrTsCap = 1000000;

results = []; % columns: L, ya, worstS11dB_in_band, minS11dB, minS11freqMHz, converged(1/0)
freq = linspace(f_low*0.7, f_high, 601);
band_idx = freq >= f_low & freq <= f_high;

for k = 1:size(combos,1)
    taper.L  = combos(k,1);
    taper.ya = combos(k,2);
    board.length = feed.throat_run + taper.L + 20; % 20 mm margin
    board.width  = 2*(taper.ya + 20);               % 20 mm margin each side

    fprintf('\n=== Combo %d/%d: taper.L=%d mm, taper.ya=%d mm (board %dx%d mm) ===\n', ...
        k, size(combos,1), taper.L, taper.ya, board.width, board.length);

    FDTD = InitFDTD('NrTs', NrTsCap);
    FDTD = SetGaussExcite(FDTD, f0, fc);
    BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'};
    FDTD = SetBoundaryCond(FDTD, BC);

    x_pts = linspace(0, taper.L, N);
    RL = taper.R * taper.L;
    C1 = (taper.ya - taper.y0) / (exp(RL) - 1);
    C2 = (taper.y0*exp(RL) - taper.ya) / (exp(RL) - 1);
    y_pts = C1*exp(taper.R*x_pts) + C2;
    y_taper_abs = feed.throat_run + x_pts;

    inner_x = y_pts; inner_y = y_taper_abs;
    outer_x = ones(1,N) * (board.width/2); outer_y = fliplr(y_taper_abs);
    Ppos = [ [inner_x, outer_x] ; [inner_y, outer_y] ];
    Pneg = [ -[inner_x, outer_x] ; [inner_y, outer_y] ];

    CSX = InitCSX();
    SimBox = [board.width+300, board.length+300, 300];
    mesh = struct();
    mesh.x = [-SimBox(1)/2 SimBox(1)/2];
    mesh.y = [-50 SimBox(2)-50];
    mesh.z = [-SimBox(3)/2 SimBox(3)/2];

    CSX = AddMaterial(CSX, 'substrate');
    CSX = SetMaterialProperty(CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa);
    start = [-board.width/2, 0, 0];
    stop  = [ board.width/2, board.length, substrate.thickness];
    CSX = AddBox(CSX, 'substrate', 0, start, stop);

    CSX = AddMetal(CSX, 'copper');
    start = [-cpw.center_w/2, 0, substrate.thickness];
    stop  = [ cpw.center_w/2, feed.throat_run, substrate.thickness];
    CSX = AddBox(CSX, 'copper', 10, start, stop);

    CSX = AddPolygon(CSX, 'copper', 10, 2, substrate.thickness, Ppos, 'CoordSystem', 0);
    CSX = AddPolygon(CSX, 'copper', 10, 2, substrate.thickness, Pneg, 'CoordSystem', 0);

    Pfeed_pos = [ [cpw.center_w/2+cpw.gap, board.width/2, board.width/2, cpw.center_w/2+cpw.gap] ; ...
                  [0, 0, feed.throat_run, feed.throat_run] ];
    Pfeed_neg = [ -[cpw.center_w/2+cpw.gap, board.width/2, board.width/2, cpw.center_w/2+cpw.gap] ; ...
                  [0, 0, feed.throat_run, feed.throat_run] ];
    CSX = AddPolygon(CSX, 'copper', 10, 2, substrate.thickness, Pfeed_pos, 'CoordSystem', 0);
    CSX = AddPolygon(CSX, 'copper', 10, 2, substrate.thickness, Pfeed_neg, 'CoordSystem', 0);

    start = [ cpw.center_w/2, 0, substrate.thickness];
    stop  = [ cpw.center_w/2+cpw.gap, 0.1, substrate.thickness];
    [CSX, port] = AddLumpedPort(CSX, 15, 1, feed.R_ref, start, stop, [1 0 0], true);

    mesh = DetectEdges(CSX, mesh, 'ExcludeProperty', {'copper'});
    mesh = DetectEdges(CSX, mesh, 'SetProperty', {'copper'}, '2D_Metal_Edge_Res', c0/(f_high)/unit/20);
    pre_tol = 0.15;
    mesh.x = sort(mesh.x); mesh.x = mesh.x([true diff(mesh.x) > pre_tol]);
    mesh.y = sort(mesh.y); mesh.y = mesh.y([true diff(mesh.y) > pre_tol]);
    mesh.z = sort(mesh.z); mesh.z = mesh.z([true diff(mesh.z) > pre_tol]);
    mesh = SmoothMesh(mesh, c0/(f_high)/unit/10);
    merge_tol = 1e-2;
    mesh.x = sort(mesh.x); mesh.x = mesh.x([true diff(mesh.x) > merge_tol]);
    mesh.y = sort(mesh.y); mesh.y = mesh.y([true diff(mesh.y) > merge_tol]);
    mesh.z = sort(mesh.z); mesh.z = mesh.z([true diff(mesh.z) > merge_tol]);
    CSX = DefineRectGrid(CSX, unit, mesh);

    start = [mesh.x(2)     mesh.y(2)     mesh.z(2)];
    stop  = [mesh.x(end-3) mesh.y(end-3) mesh.z(end-3)];
    [CSX, nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);

    Sim_Path = sprintf('Sweep_L%d_ya%d', taper.L, taper.ya);
    Sim_CSX  = 'VIV.xml';
    [status, message, messageid] = rmdir(Sim_Path, 's');
    [status, message, messageid] = mkdir(Sim_Path);
    WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);

    try
        RunOpenEMS(Sim_Path, Sim_CSX, '--debug-PEC');
    catch err
        fprintf('  -> run failed/interrupted: %s\n', err.message);
    end

    try
        port = calcPort(port, Sim_Path, freq);
        s11dB = 20*log10(abs(port.uf.ref ./ port.uf.inc));
        worst = max(s11dB(band_idx));
        [best_val, best_ix] = min(s11dB);
        converged = 1; % we don't parse the solver log here; verify manually for the winner
        results = [results; taper.L taper.ya worst best_val freq(best_ix)/1e6 converged]; %#ok<AGROW>
        fprintf('  -> worst-case S11 in 850-2000 MHz band: %.2f dB | deepest dip: %.2f dB @ %.0f MHz\n', ...
            worst, best_val, freq(best_ix)/1e6);
    catch err
        fprintf('  -> post-processing failed: %s\n', err.message);
        results = [results; taper.L taper.ya NaN NaN NaN 0]; %#ok<AGROW>
    end
end

fprintf('\n\n========== SWEEP SUMMARY (sorted best-first) ==========\n');
fprintf('%8s %8s %18s %14s %12s\n', 'taper.L', 'taper.ya', 'worst S11 (dB)', 'deepest dip', 'at freq(MHz)');
[~, order] = sort(results(:,3)); % ascending worst-case S11 (most negative = best)
for i = 1:size(results,1)
    r = results(order(i),:);
    fprintf('%8d %8d %18.2f %14.2f %12.0f\n', r(1), r(2), r(3), r(4), r(5));
end

best = results(order(1),:);
fprintf('\nBEST combination: taper.L=%d mm, taper.ya=%d mm (worst-case S11 = %.2f dB across 850-2000 MHz)\n', ...
    best(1), best(2), best(3));
fprintf('IMPORTANT: this used the coarse exploration mesh and a reduced NrTs cap (%d) for speed.\n', NrTsCap);
fprintf('Re-run this exact combination in Vivaldi_OpenEMS_Start.m (fine mesh, full NrTs) to confirm\n');
fprintf('before treating it as your final design.\n');
