% wg_alumina_slotted_gain.m
% Alumina-filled slotted waveguide (openEMS) with speed/final toggle
% + S11 + Far-field + Directivity + (Realized) Gain.
% Axes: X=b (thickness), Y=a (broad), Z=length (propagation).

clear; clc;

%% ---------------- Toggles ----------------
speedMode      = true;     % set false for full/final run
forcePatterns  = true;     % set true to compute patterns even in FAST mode

%% ---------------- User Parameters (mm) ----------------
f0      = 10.5e9;                   % Hz
fbw     = 0.20;                     % S11 band ±10%
eps_r   = 9.8;                      % alumina
tanD    = 1e-4;                     % alumina loss tangent

a       = 8.5;                      % Y (mm)
b       = 3.5;                      % X (mm)
twall   = 0.6;                      % wall thickness (mm)

% slot rules
slotWidth_rule = 0.02;              % w ~ 0.02*lambda0
slotLen_rule   = 0.47;              % L ~ 0.47*lambda_g
x_off_baseFrac = 0.25;              % base offset ~ 0.25*a
x_off_minFrac  = 0.22;              % clamp min
x_off_maxFrac  = 0.30;              % clamp max
kaiser_beta    = 1.65;              % ~ -25 dB-like taper

% coax (SMA-ish) defaults
r_in_default  = 0.30;               % mm
r_out_default = 1.50;               % mm

% speed/final presets
if speedMode
    Nslots = 9;
    r_in   = 0.50;                  % thicker pin -> coarser mesh
    r_out  = 1.80;
    slot_W_override = 0.8;          % widen slot to 0.8 mm
    endCrit = 5e-4;
    ratio   = 2.2;
    doNF2FF = false;
else
    Nslots = 32;
    r_in   = r_in_default;
    r_out  = r_out_default;
    slot_W_override = [];
    endCrit = 1e-4;
    ratio   = 1.7;
    doNF2FF = true;
end
% allow forcing patterns on in fast mode
if forcePatterns
    doNF2FF = true;
end

r_os = r_out + 0.40;                % outer shell radius (mm)

%% ---------------- Derived guide quantities ----------------
lambda0 = (3e8/f0)*1e3;             % mm
fc10    = (3e8)/(2*(a*1e-3)*sqrt(eps_r));
assert(f0>fc10, 'f0 below TE10 cutoff with alumina fill.');

lambda_g = lambda0 / sqrt(1 - (fc10/f0)^2);
dz       = lambda_g/2;
slot_L   = slotLen_rule   * lambda_g;
slot_W   = slotWidth_rule * lambda0;
if ~isempty(slot_W_override), slot_W = slot_W_override; end
x_off_base = x_off_baseFrac * a;

% longitudinal layout (Z)
z_first  = -((Nslots-1)/2)*dz;
z_feed   =  z_first - lambda_g/4;
z_short  =  z_feed  - lambda_g/4;

% margins along Z
if speedMode
    front_margin = dz;
    back_margin  = dz;
else
    front_margin = 2*dz;
    back_margin  = 2*dz;
end

z_min = z_short - back_margin;
z_max = -z_first + front_margin;
Lz    = (z_max - z_min);

fprintf('lambda0=%.2f mm, lambda_g=%.2f mm, dz=%.2f mm\n',lambda0,lambda_g,dz);
fprintf('z_first=%.2f mm, z_feed=%.2f mm, z_short=%.2f mm\n',z_first,z_feed,z_short);

%% ---------------- openEMS / CSX ----------------
unit = 1e-3;                        % mm
FDTD = InitFDTD('EndCriteria', endCrit);
FDTD = SetGaussExcite(FDTD, f0, (f0*fbw)/2);
FDTD = SetBoundaryCond(FDTD, {'PML_8','PML_8','PML_8','PML_8','PML_8','PML_8'});

CSX = InitCSX();

% materials
CSX = AddMetal(CSX,'PEC');

CSX = AddMaterial(CSX,'ALUMINA');
CSX = SetMaterialProperty(CSX,'ALUMINA','Epsilon',eps_r,'Mue',1, ...
    'Kappa', 2*pi*f0*8.854187817e-12*eps_r*tanD);

CSX = AddMaterial(CSX,'AIR');
CSX = SetMaterialProperty(CSX,'AIR','Epsilon',1,'Mue',1);

CSX = AddMaterial(CSX,'PTFE');
CSX = SetMaterialProperty(CSX,'PTFE','Epsilon',2.1);

%% ---------------- Geometry (X=b, Y=a, Z=length) ----------------
x0 = -b/2;  x1 =  b/2;
y0 = -a/2;  y1 =  a/2;
z0 = z_min; z1 = z_max;

% alumina fill
CSX = AddBox(CSX,'ALUMINA',1,[x0 y0 z0],[x1 y1 z1]);

% PEC walls + backshort
CSX = AddBox(CSX,'PEC',10,[x1           y0 z0],[x1+twall  y1 z1]); % top
CSX = AddBox(CSX,'PEC',10,[x0-twall     y0 z0],[x0        y1 z1]); % bottom
CSX = AddBox(CSX,'PEC',10,[x0 y1        z0],[x1 y1+twall  z1]);    % +Y
CSX = AddBox(CSX,'PEC',10,[x0 y0-twall  z0],[x1 y0        z1]);    % -Y
CSX = AddBox(CSX,'PEC',12,[x0 y0 z_short],[x1 y1 z_short+twall]);  % backshort

% Kaiser taper -> slot offsets (alternate signs for dz=lambda_g/2)
w = kaiser(Nslots, kaiser_beta)';  w = w/max(w);
x_off_mag = x_off_base * sqrt(w);
x_off_mag = min(max(x_off_mag, x_off_minFrac*a), x_off_maxFrac*a);
signs     = (-1).^(0:Nslots-1);
y_centers = signs .* x_off_mag;

% Drill each slot through top wall (AIR with higher priority)
for n = 1:Nslots
    zc  = z_first + (n-1)*dz;
    y_c = y_centers(n);
    y1s = y_c - slot_W/2;
    y2s = y_c + slot_W/2;
    z1s = zc  - slot_L/2;
    z2s = zc  + slot_L/2;
    CSX = AddBox(CSX,'AIR',20,[x1 y1s z1s],[x1+twall y2s z2s]);
end

%% ---------------- Mesh (STRUCT; BEFORE AddCoaxialPort) ----------------
if speedMode
    maxCell = min([slot_W, r_in, lambda0/15]);
else
    maxCell = min([slot_W, r_in, lambda0/20]);
end

% Y lines (include slot edges)
dy_lines = unique(sort([y0-2, y0, y1, y1+2, y_centers-slot_W/2, y_centers+slot_W/2]));
dy = SmoothMeshLines(dy_lines, maxCell, ratio);

% Z lines (include slot edges + feed/short)
slot_z_edges = zeros(1, 2*Nslots);
for n = 1:Nslots
    zc = z_first + (n-1)*dz;
    slot_z_edges(2*n-1) = zc - slot_L/2;
    slot_z_edges(2*n)   = zc + slot_L/2;
end
dz_lines = unique(sort([z0-5, z0, slot_z_edges, z_feed, z_short, z1, z1+5]));
dz = SmoothMeshLines(dz_lines, maxCell, ratio);

% X lines (thickness)
dx = SmoothMeshLines([x0-2, x0, x1, x1+2], maxCell, ratio);

mesh = struct('x', dx, 'y', dy, 'z', dz);
CSX  = DefineRectGrid(CSX, unit, mesh);

%% ---------------- Coaxial probe feed (AFTER mesh) ----------------
% Coax along +Y from -Y wall at z = z_feed
start = [0, y0 - twall - 0.5, z_feed];   % outside PEC by 0.5 mm
stop  = [0, y0 + 0.05,        z_feed];   % 50 µm inside alumina

% 11-arg API: AddCoaxialPort(CSX,prio,portnr,'PEC','PTFE',start,stop,dir,r_i,r_o,r_os,...)
[CSX, port] = AddCoaxialPort( ...
    CSX, 100, 1, 'PEC', 'PTFE', start, stop, 'y', r_in, r_out, r_os, ...
    'ExciteAmp', 1.0);

%% ---------------- NF2FF box (robust across versions) ----------------
doNF2FF = doNF2FF || forcePatterns;   % keep your toggles

if doNF2FF
    nf_box_name = 'nf2ff';
    nf_start = [x0 - (twall+1), y0 - (twall+1), z0 - 5];
    nf_stop  = [x1 + (twall+1), y1 + (twall+1), z1 + 5];

    % Try modern creator that RETURNS a struct:
    nf2ff = [];
    try
        nf2ff = CreateNF2FFBox(CSX, nf_box_name, nf_start, nf_stop);
    catch
        % Fallback: older helper that RETURNS nothing; synthesize the struct
        try
            AddNF2FFBox(CSX, nf_box_name, nf_start, nf_stop);
            nf2ff = struct('name', nf_box_name);
        catch
            % Some very old versions add via AddBox/SetNF2FFBox—last resort:
            error('NF2FF box creation failed. Your openEMS build lacks Create/AddNF2FFBox.');
        end
    end
end



%% ---------------- Write / View / Run ----------------
if speedMode, Sim_Path = 'sim_wg_alumina_fast';
else,        Sim_Path = 'sim_wg_alumina_final'; end
Sim_CSX  = 'wg_alumina.xml';
[~,~] = rmdir(Sim_Path, 's'); mkdir(Sim_Path);
WriteOpenEMS(fullfile(Sim_Path, Sim_CSX), FDTD, CSX);

if exist('AppCSXCAD','file'), AppCSXCAD(fullfile(Sim_Path, Sim_CSX)); end
RunOpenEMS(Sim_Path, Sim_CSX, '--engine=multithreaded');

%% ---------------- S11 (robust block) ----------------
f1 = f0*(1 - fbw/2);  
f2 = f0*(1 + fbw/2);
Npts = 401;
fvec = linspace(f1, f2, Npts);

try
    port = calcPort(port, Sim_Path, fvec, 'RefImpedance', 50);
catch
    port = calcPort(port, Sim_Path, fvec);
end

if isfield(port, 'uf') && isfield(port.uf,'ref') && isfield(port.uf,'inc')
    S11 = port.uf.ref ./ port.uf.inc;
elseif isfield(port, 'U') && isfield(port.U,'ref') && isfield(port.U,'inc')
    S11 = port.U.ref ./ port.U.inc;
else
    error('Unknown port struct fields. Inspect "port".');
end

figure; plot(fvec/1e9, 20*log10(abs(S11)), 'LineWidth',1.5);
grid on; xlabel('Frequency (GHz)'); ylabel('|S_{11}| (dB)');
if speedMode, title(sprintf('|S_{11}| alumina %d-slot (FAST)', Nslots));
else,        title(sprintf('|S_{11}| alumina %d-slot (FINAL)', Nslots)); end

% S11 at f0 for mismatch efficiency
[~, idx0] = min(abs(fvec - f0));
Gamma0 = S11(idx0);
eta_m  = 1 - abs(Gamma0)^2;

%% ---------------- Far-field, Directivity, Gain ----------------
if doNF2FF
    theta = linspace(0,180,361);   % deg
    phi   = linspace(0,360,721);   % deg

    % Ensure we pass a struct with .name; if a string slipped through, wrap it
    if ischar(nf2ff)
        nf2ff = struct('name', nf2ff);
    end

    try
        ff3d = CalcNF2FF(nf2ff, Sim_Path, f0, theta, phi, 'OutDir','nf2ff');
    catch
        % Some builds don’t support 'OutDir' option
        ff3d = CalcNF2FF(nf2ff, Sim_Path, f0, theta, phi);
    end


    % Radiation intensity U = |E|^2/(2*eta0) at r=1 m
    eta0 = 376.730313668;
    U = (E.^2) ./ (2*eta0);      % W/sr

    % Integrate over sphere to get Prad
    th = deg2rad(theta(:));      % Nx1
    ph = deg2rad(phi(:)).';      % 1xM
    dth = [diff(th); th(end)-th(end-1)];
    dph = [diff(ph), ph(end)-ph(end-1)];
    [DTH, DPH] = ndgrid(dth, dph);
    [SINE, ~]  = ndgrid(sin(th), ones(size(ph)));
    dOmega = SINE .* DTH .* DPH;
    Prad = sum(U .* dOmega, 'all');

    % Directivity
    D = 4*pi * U / Prad;
    [Dmax, imax] = max(D(:));
    [i_th, i_ph] = ind2sub(size(D), imax);
    th_max = theta(i_th);  ph_max = phi(i_ph);

    % Accepted power Pacc from port struct (best-effort across versions)
    Pacc = NaN;
    if isfield(port,'P_inc') && isfield(port,'P_ref')
        Pinc0 = port.P_inc(idx0); Pref0 = port.P_ref(idx0);
        Pacc = real(Pinc0 - Pref0);
    elseif isfield(port,'uf') && isfield(port,'if') ...
           && isfield(port.uf,'inc') && isfield(port.if,'inc') ...
           && isfield(port.uf,'ref') && isfield(port.if,'ref')
        Uinc = port.uf.inc(idx0);  Iinc = port.if.inc(idx0);
        Uref = port.uf.ref(idx0);  Iref = port.if.ref(idx0);
        Pacc = real(0.5*(Uinc*conj(Iinc) - Uref*conj(Iref)));
    elseif isfield(port,'U') && isfield(port,'I') ...
           && isfield(port.U,'inc') && isfield(port.I,'inc') ...
           && isfield(port.U,'ref') && isfield(port.I,'ref')
        Uinc = port.U.inc(idx0);  Iinc = port.I.inc(idx0);
        Uref = port.U.ref(idx0);  Iref = port.I.ref(idx0);
        Pacc = real(0.5*(Uinc*conj(Iinc) - Uref*conj(Iref)));
    end
    if ~isfinite(Pacc) || Pacc<=0
        warning('Pacc not found in port struct; assuming Pacc = 1 W.');
        Pacc = 1.0;
    end
    eta_rad = Prad / Pacc;  eta_rad = max(min(eta_rad,1),0);

    % Gains
    Gmax      = eta_rad * Dmax;
    Greal_max = eta_m  * Gmax;
    Gmax_dBi      = 10*log10(Gmax);
    Greal_max_dBi = 10*log10(Greal_max);

    fprintf('\n=== Gain summary @ %.3f GHz ===\n', f0/1e9);
    fprintf('Prad = %.3f W, Pacc = %.3f W  =>  eta_rad = %.1f %%\n', Prad, Pacc, 100*eta_rad);
    fprintf('Mismatch efficiency eta_m = %.1f %%  (|S11|=%.2f dB)\n', 100*eta_m, 20*log10(abs(Gamma0)));
    fprintf('Dmax = %.2f (%.2f dBi) at (theta=%.1f°, phi=%.1f°)\n', Dmax, 10*log10(Dmax), th_max, ph_max);
    fprintf('Gmax      = %.2f dBi\n', Gmax_dBi);
    fprintf('Greal_max = %.2f dBi (includes mismatch)\n\n', Greal_max_dBi);

    % Quick E/H-plane plots (normalized)
    [~,i_phi0] = min(abs(phi-0));
    EnormE = E(:,i_phi0) / max(E(:,i_phi0));
    figure; plot(theta, 20*log10(abs(EnormE)+1e-15),'LineWidth',1.5);
    grid on; xlabel('\theta (deg)'); ylabel('Normalized |E| (dB)'); ylim([-60 0]);
    title('E-plane cut (\phi \approx 0^\circ)');

    [~,i_phiH] = min(abs(phi-90));
    EnormH = E(:,i_phiH) / max(E(:,i_phiH));
    figure; plot(theta, 20*log10(abs(EnormH)+1e-15),'LineWidth',1.5);
    grid on; xlabel('\theta (deg)'); ylabel('Normalized |E| (dB)'); ylim([-60 0]);
    title('H-plane cut (\phi \approx 90^\circ)');
end

