% Vivaldi_OpenEMS_Start.m
%
% Starter model for a CPW-fed, single-copper-layer, planar broadband
% Vivaldi antenna on FR4, for the Studienprojekt "Design einer
% Breitband-Vivaldi-Antenne" (Prof. Passoke).
%
% This script reuses the overall structure of the supplied Dipol_Passoke.m
% (setup / substrate / metal / port / mesh / NF2FF / run / evaluate) and
% replaces the straight dipole arms with an exponentially tapered Vivaldi
% flare, built with AddPolygon (the P2 example that was commented out in
% Dipol_Passoke.m is the template this is based on).
%
% STATUS: first-guess starting geometry, NOT yet optimized. Dimensions
% come from the "Vivaldi antenna initial design" section of the project
% strategy document. You must run this, check S11/Zin, and then sweep
% parameters (Section 6 of the strategy doc) to converge on a design with
% S11 < -10 dB from 850 MHz upward.
%
% ASSUMPTIONS made here that you should confirm with Prof. Passoke:
%   - Design/simulation band assumed 850 MHz - 3 GHz (upper limit not
%     given in the task sheet).
%   - Max board size assumed to fit within ~200 mm x 150 mm.
%   - CPW-fed single-layer design chosen instead of a double-sided
%     microstrip-to-slotline balun (task says a balun "may" be needed,
%     not that it is mandatory) - confirm this choice is acceptable.

close all
clear
clc

%% ------------------------------------------------------------------
%% 1) Constants & frequency setup
%% ------------------------------------------------------------------
physical_constants;
unit = 1e-3; % all lengths below are in mm

f_low  = 850e6;   % required lower band edge
f_high = 2000e6;  % ASSUMED upper band edge -> confirm with professor

f0 = (f_low + f_high)/2;   % Gaussian excitation center frequency
fc = (f_high - f_low)/2 + 0.3e9; % 20 dB corner freq, with margin below f_low
                                  % so the excitation spectrum still has
                                  % useful energy right down at 850 MHz

FDTD = InitFDTD('NrTs', 3500000); % INCREASED from 2000000: the last run
                                   % (taper.L=160, ya=70, R=0.014) hit this
                                   % cap before reaching -50dB decay.
                                   % timesteps than a simple dipole to
                                   % decay fully - watch the end criteria
                                   % in the solver log; increase further
                                   % if it stops before decaying enough
FDTD = SetGaussExcite(FDTD, f0, fc);

BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % start with MUR like the
                                             % dipole example; switch to
                                             % PML_8 if you see far-field
                                             % artifacts from reflections
FDTD = SetBoundaryCond(FDTD, BC);

%% ------------------------------------------------------------------
%% 2) Substrate parameters
%% ------------------------------------------------------------------
% Air (kappa computed at design center frequency, same pattern as the
% dipole starter file)
luft.epsR  = 1;
luft.kappa = 0; % lossless air

% FR4 substrate: epsR and tan delta from the datasheet / task sheet.
% kappa = tan(delta) * 2*pi*f*eps0*epsR  (constant-kappa approximation,
% evaluated at the band center f0 - valid over a moderate bandwidth;
% mention this approximation explicitly in your report)
substrate.epsR      = 4.3;
substrate.tand      = 0.023;
substrate.kappa     = substrate.tand * 2*pi*f0 * substrate.epsR * EPS0;
substrate.thickness = 1;      % mm, given by the task/datasheet

%% ------------------------------------------------------------------
%% 3) Vivaldi geometry parameters (first-guess values, see strategy doc)
%% ------------------------------------------------------------------
board.length = 275;   % mm, taper direction (y-axis here) -- fits taper.L=240
                       % (throat_run 15 + 240 = 255, +20 mm margin).
board.width  = 160;   % mm, transverse direction (x-axis here) -- fits
                       % taper.ya=60 (2*(60+20) mm margin).

feed.throat_run = 15;   % mm, straight CPW section before the taper starts
feed.R_ref      = 50;   % ohm, target system/port impedance (SMA/CPW)

% CPW dimensions defined BEFORE the taper block on purpose: taper.y0 below
% is derived from these so the ground plane's inner edge is geometrically
% continuous at the throat/taper junction (see bug note below).
cpw.center_w = 3.0;  % mm, CPW center conductor width (feed section)
cpw.gap      = 0.4;  % mm, CPW gap, center conductor to ground
% NOTE: cpw.center_w / cpw.gap are analytic-formula starting guesses for
% ~50 ohm on 1 mm FR4 (er=4.3). Finite ground width and the taper make
% analytic CPW formulas inaccurate here - verify/tune these two numbers
% against the simulated Zin, don't trust them blindly.

taper.L  = 240;   % mm, taper length (throat -> aperture) -- BEST result
                  % from the automated sweep (Vivaldi_Sweep.m): 6 combos
                  % of taper.L in {160,200,240} x taper.ya in {50,60} were
                  % tested. L=240/ya=60 had the best (least bad) worst-
                  % case S11 across 850-2000 MHz (-3.18 dB), though NONE
                  % of the 6 combos reached the -10 dB target continuously.
                  % This is the production run to confirm that sweep
                  % result with a finer mesh and full NrTs.
% BUG FIX (see Report_Figures / KI-Nutzung notes): taper.y0 was previously
% hardcoded to 1.0 mm, independent of cpw.center_w/cpw.gap. Since the feed
% section's ground plane starts at x = cpw.center_w/2 + cpw.gap = 1.9 mm,
% but the taper's ground inner edge started at x = 1.0 mm at the very same
% y-position (y = feed.throat_run), the ground copper jumped inward past
% the end of the 1.5 mm-half-width center conductor at that seam - i.e.
% the ground plane physically overlapped/shorted the center conductor's
% open end right at the throat/taper junction. This is almost certainly
% the root cause of the S11 staying far from -10 dB and of the simulation
% instabilities (oscillating/growing energy) seen when other dimensions
% were changed, since changing board.length/taper.L moves the mesh lines
% around this short and changes how it resonates.
% FIX: start the taper's slot exactly where the feed section's gap ends,
% so the ground plane is continuous (no overlap, no unintended gap).
taper.y0 = cpw.center_w/2 + cpw.gap;  % mm, HALF throat gap width at x=0
                                       % (start of taper) = 1.9 mm
taper.ya = 60.0;  % mm, HALF aperture width at x=L (end of taper) --
                  % see taper.L note above (best sweep result).
taper.R  = 0.023; % 1/mm, taper rate - KEY OPTIMIZATION PARAMETER --
                  % REVERTED to the confirmed-stable value. Not touching
                  % this again; taper.L is the lever being pushed now.

%% ------------------------------------------------------------------
%% 4) Metal polygons: exponential taper curve generation
%% ------------------------------------------------------------------
% y(x) = C1*exp(R*x) + C2 ,  0 <= x <= L
% with y(0) = y0 and y(L) = ya (see strategy doc Section 5 for derivation)
N = 10; % kept at 10 (proven stable at L=240/ya=60 in the sweep) - not
        % risking a SmoothMesh hang on the production run by raising this.
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;

% Place the throat at y-position 'feed.throat_run' measured from the
% board's feed edge (y=0), so the full antenna occupies
% y = 0 ... feed.throat_run + taper.L along the board.
y_taper_abs = feed.throat_run + x_pts; % absolute y-coordinate along board

% Build the closed polygon for ONE ground wing (the +x side): the ground
% plane fills the board from the taper edge out to the board's side edge,
% from the throat all the way to the aperture end. Point lists are
% assembled as [x-row ; y-row], matching the AddPolygon(...,'CoordSystem',0)
% convention used in the commented-out P2 example of Dipol_Passoke.m.

% Inner boundary (the slot edge), going from throat to aperture:
inner_x =  y_pts;                 % x-coordinate = taper half-width
inner_y =  y_taper_abs;            % y-coordinate = position along board

% Outer boundary (board edge), going back from aperture to throat:
outer_x =  ones(1,N) * (board.width/2);
outer_y =  fliplr(y_taper_abs);

Ppos = [ [inner_x, outer_x] ; [inner_y, outer_y] ]; % ground wing, +x side
Pneg = [ -[inner_x, outer_x] ; [inner_y, outer_y] ]; % ground wing, -x side (mirrored)

%% ------------------------------------------------------------------
%% 5) CSXCAD geometry
%% ------------------------------------------------------------------
CSX = InitCSX();

% mesh air-box: must enclose the (larger) Vivaldi aperture + radiating
% space; recompute/enlarge if you change board size
SimBox = [board.width+300, board.length+300, 300];

mesh.x = [-SimBox(1)/2 SimBox(1)/2];
mesh.y = [-50 SimBox(2)-50];
mesh.z = [-SimBox(3)/2 SimBox(3)/2];

% substrate
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);

% copper: CPW feed line (straight section from y=0 to throat)
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);

% copper: ground wings (tapered), both sides, built from the polygons above
CSX = AddPolygon(CSX, 'copper', 10, 2, substrate.thickness, Ppos, 'CoordSystem', 0);
CSX = AddPolygon(CSX, 'copper', 10, 2, substrate.thickness, Pneg, 'CoordSystem', 0);

% also fill the ground on either side of the straight feed section
% (from the CPW gap out to the board edge), so the CPW line has ground
% right from the feed edge, not just starting at the throat
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);

%% ------------------------------------------------------------------
%% 6) Feed / port
%% ------------------------------------------------------------------
% Lumped port bridging the CPW gap at the feed edge (y=0), same mechanism
% as the gap-bridging port in Dipol_Passoke.m, just oriented across the
% CPW gap instead of along the dipole arms.
start = [ cpw.center_w/2, 0, substrate.thickness];
stop  = [ cpw.center_w/2+cpw.gap, 0.1, substrate.thickness]; % bridging the CPW gap
[CSX, port] = AddLumpedPort(CSX, 15, 1, feed.R_ref, start, stop, [1 0 0], true);

%% ------------------------------------------------------------------
%% 6b) Geometry sanity check (native Octave plot - no AppCSXCAD needed)
%% ------------------------------------------------------------------
% AppCSXCAD/QCSXCAD is not installed on this Mac, so CSXGeomPlot() cannot
% be used. This 2D plot draws the same copper polygons directly from the
% arrays already computed above (feed line box, ground taper wings, feed
% ground fill, and the port location) so you can catch overlap/gap bugs
% BEFORE spending 10-20+ minutes on a full FDTD run. ALWAYS look at this
% plot after changing any geometry parameter.
figure('Name', 'Geometry check (copper layer, top view)');
hold on
% feed line center conductor
fill([-cpw.center_w/2 cpw.center_w/2 cpw.center_w/2 -cpw.center_w/2], ...
     [0 0 feed.throat_run feed.throat_run], [0.85 0.5 0.2]);
% ground: feed-section fill + taper wings, both sides
fill(Pfeed_pos(1,:), Pfeed_pos(2,:), [0.7 0.7 0.7]);
fill(Pfeed_neg(1,:), Pfeed_neg(2,:), [0.7 0.7 0.7]);
fill(Ppos(1,:), Ppos(2,:), [0.7 0.7 0.7]);
fill(Pneg(1,:), Pneg(2,:), [0.7 0.7 0.7]);
% port location (bridges the CPW gap at y=0)
plot([cpw.center_w/2 cpw.center_w/2+cpw.gap], [0.05 0.05], 'r-', 'LineWidth', 4);
% board outline
plot([-board.width/2 board.width/2 board.width/2 -board.width/2 -board.width/2], ...
     [0 0 board.length board.length 0], 'k--');
axis equal
xlabel('x / mm'); ylabel('y / mm');
title('Geometry check: orange=feed conductor, grey=ground, red=port - look for overlaps/gaps');
grid on
drawnow
disp('>>> CHECK THE GEOMETRY FIGURE NOW. Only continue if the orange feed');
disp('    line and the grey ground planes do NOT overlap, and the red port');
disp('    bar sits cleanly inside the gap between them.');
% NOTE: this check now runs BEFORE mesh generation (moved here on purpose)
% so you see it quickly, without waiting through SmoothMesh, which can be
% slow/hang for larger geometries.

%% ------------------------------------------------------------------
%% 7) Mesh
%% ------------------------------------------------------------------
% EXPLORATION-PHASE MESH: deliberately coarser than a final production
% mesh (divisors reduced from /50,/20 to /20,/10) per the project's own
% Schnellplan ("use a coarser mesh for exploratory sweep runs"). Larger
% board/aperture sizes were taking 25+ minutes without fully converging
% at the finer resolution - this trades some accuracy for much faster
% iteration while we search for the right taper.L/ya/R combination.
% IMPORTANT: once dimensions are chosen, switch these back to /50 and /20
% (or finer) and re-run once for the final, report-quality result.
mesh = DetectEdges(CSX, mesh, 'ExcludeProperty', {'copper'});
mesh = DetectEdges(CSX, mesh, 'SetProperty', {'copper'}, '2D_Metal_Edge_Res', c0/(f_high)/unit/35);
% resolution raised from /20 to /35 for this production/confirmation run
% (finer than the exploration sweep, still short of the original /50 to
% avoid re-triggering a SmoothMesh hang at this larger board size).

% ROBUSTNESS FIX: pre-thin the raw DetectEdges line set BEFORE calling
% SmoothMesh. SmoothMesh's recursive line-smoothing search has repeatedly
% hung (100% CPU, 10-30+ min, zero output) whenever the taper curve's N
% sample points happen to land close together in absolute position -
% this has now happened with several different, unrelated parameter
% changes (taper.L, board.width, taper.R), so it is SmoothMesh's own
% algorithm that is fragile against dense/close line clusters, not any
% one specific bug in our geometry. Thinning the raw line set first (any
% two lines closer than pre_tol mm are merged) removes those clusters
% before the expensive recursive step ever sees them.
pre_tol = 0.15; % mm (still well below the smallest real feature, the 0.4 mm CPW gap)
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/15); % raised from /10 to /15

% BUG FIX: floating-point rounding in the exponential taper formula can
% produce two mesh lines that are meant to be identical (e.g. both at
% x = cpw.center_w/2+cpw.gap) but differ by ~1e-16 mm due to roundoff.
% openEMS then builds a real (if absurdly tiny) cell between them, which
% collapses the FDTD timestep to ~1e-20 s and makes the run meaningless.
% Fix: snap/merge any mesh lines closer than 1e-6 mm (still far finer than
% any real feature) on each axis before handing the mesh to openEMS.
merge_tol = 1e-2; % mm (10 micron; still 40x below the smallest real feature, the 0.4 mm gap)
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);

%% ------------------------------------------------------------------
%% 8) NF2FF box
%% ------------------------------------------------------------------
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);

%% ------------------------------------------------------------------
%% 9) Run simulation
%% ------------------------------------------------------------------
Sim_Path = 'Vivaldi01';
Sim_CSX  = 'VIV01.xml';

[status, message, messageid] = rmdir(Sim_Path, 's');
[status, message, messageid] = mkdir(Sim_Path);

WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);
% CSXGeomPlot([Sim_Path '/' Sim_CSX]); % ALWAYS check this plot before running -
                                      % confirm the taper/ground/feed look
                                      % correct before spending simulation time

openEMS_opts = '--debug-PEC';
RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts);

%% ------------------------------------------------------------------
%% 10) Post-processing: S11 / Zin
%% ------------------------------------------------------------------
freq = linspace(f_low*0.7, f_high, 601); % start below f_low to see the
                                          % full roll-off of the match
port = calcPort(port, Sim_Path, freq);

Zin = port.uf.tot ./ port.if.tot;
s11 = port.uf.ref ./ port.uf.inc;

figure
plot(freq/1e6, real(Zin), 'k-', 'Linewidth', 2); hold on; grid on
plot(freq/1e6, imag(Zin), 'r--', 'Linewidth', 2);
title('Feed point impedance');
xlabel('frequency f / MHz'); ylabel('impedance Z_{in} / Ohm');
legend('real', 'imag');

figure
plot(freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2); grid on
hold on
plot([f_low f_low]/1e6, ylim, 'b--'); % marker at the required 850 MHz edge
yline(-10, 'g--'); % -10 dB target line
title('Reflection coefficient S_{11}');
xlabel('frequency f / MHz'); ylabel('|S_{11}| / dB');

drawnow

%% ------------------------------------------------------------------
%% 11) Post-processing: far-field (NF2FF), directivity, efficiency
%% ------------------------------------------------------------------
f_res_ind = find(20*log10(abs(s11)) == min(20*log10(abs(s11))));
f_res = freq(f_res_ind);

disp('calculating far field at phi=[0 90] deg...');
nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, [-180:2:180]*pi/180, [0 90]*pi/180);

disp(['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']);
disp(['directivity: Dmax = ' num2str(nf2ff.Dmax) ' (' num2str(10*log10(nf2ff.Dmax)) ' dBi)']);
disp(['efficiency: nu_rad = ' num2str(100*nf2ff.Prad./port.P_inc(f_res_ind)) ' %']);

figure
polarFF(nf2ff, 'xaxis', 'theta', 'param', [1 2], 'normalize', 1);

figure
plotFFdB(nf2ff, 'xaxis', 'theta', 'param', [1 2]);

drawnow

disp('calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...');
thetaRange = (0:2:180);
phiRange   = (0:2:360) - 180;
nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180, 'Verbose', 1, 'Outfile', '3D_Pattern.h5');

figure
plotFF3D(nf2ff, 'logscale', -20);

E_far_normalized = nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)) * nf2ff.Dmax;
DumpFF2VTK([Sim_Path '/3D_Pattern.vtk'], E_far_normalized, thetaRange, phiRange, 'scale', 1e-3);

%% ------------------------------------------------------------------
%% NEXT STEPS (see strategy doc Section 6 for full detail)
%% ------------------------------------------------------------------
% 1. Check CSXGeomPlot output BEFORE trusting the simulation - confirm
%    the taper shape, ground wings, and feed gap look geometrically
%    sensible (no overlaps, no gaps you didn't intend).
% 2. Look at the S11 plot: find the lowest frequency where it crosses
%    below -10 dB. If that's above 850 MHz, increase taper.L and/or
%    taper.ya first (biggest levers), then re-run.
% 3. Once S11 looks reasonable, verify Zin real part is close to 50 ohm
%    and imaginary part is small, over your target band.
% 4. Sweep taper.R, cpw.center_w, cpw.gap, feed.throat_run one at a time,
%    tracking how the -10 dB crossing frequency and Zin move.
% 5. After the design converges, refine the mesh (smaller
%    2D_Metal_Edge_Res / SmoothMesh target) and re-run once to confirm
%    the result is mesh-converged, not a meshing artifact.
