Filter design using fractral geometry (Sierpinski resonator with a microstrip feed)

% Clear and initialize

close all; clc;

% Add OpenEMS and CSXCAD paths

addpath('C:\openEMS\matlab');

addpath('C:\openEMS\CSXCAD\matlab');

% Constants

c0 = 299792458;  % Speed of light in vacuum [m/s]

unit = 1e-3;     % millimeter

%% Adjustable Parameters

iter = 2;                    % Sierpinski depth

sideLength = 30;            % Triangle side [mm]

substrate_thickness = 1.6;  % mm

eps_r = 4.4;                % Dielectric constant of FR4

f_start = 1e9;              % Simulation start frequency

f_stop = 10e9;              % Simulation stop frequency

res = 30;                   % Mesh resolution factor

%% Calculate microstrip width for 50 Ohm

Z0 = 50; h = substrate_thickness; er = eps_r;

if Z0 <= 60

    W = (8 * h) * exp(Z0 / 60 * sqrt((eps_r + 1)/2)) / ...

        (exp(2 * Z0 / 60 * sqrt((eps_r + 1)/2)) - 2);

else

    B = (377 * pi) / (2 * Z0 * sqrt(eps_r));

    W = 2 * h / pi * (B - 1 - log(2*B - 1) + ...

        ((eps_r - 1)/(2*eps_r)) * (log(B - 1) + 0.39 - 0.61/eps_r));

end

%% Simulation box and mesh

SimBox = [sideLength*1.5, sideLength*1.5, 10]; % mm

max_res = c0 / (f_stop * res) / unit;

%% Create CSX geometry

CSX = InitCSX();

CSX = AddMaterial(CSX, 'FR4');

CSX = SetMaterialProperty(CSX, 'FR4', 'Epsilon', eps_r);

CSX = AddBox(CSX, 'FR4', 0, [0 0 0], [SimBox(1) SimBox(2) h]);

% Ground plane

CSX = AddMetal(CSX, 'Ground');

CSX = AddBox(CSX, 'Ground', 0, [0 0 0], [SimBox(1) SimBox(2) 0]);

% Metal layer for microstrip + resonator

CSX = AddMetal(CSX, 'Patch');

% Microstrip feed line

feed_length = sideLength * 0.5;

feed_width = W;

feed_x = sideLength / 2 - feed_length/2;

CSX = AddBox(CSX, 'Patch', 0, ...

    [feed_x, 0, h], [feed_x + feed_length, feed_width, h]);

%% Triangle vertices (attach to feed end)

tx = feed_x + feed_length;

ty = feed_width/2;

h_tri = sqrt(3)/2 * sideLength;

v1 = [tx, ty];

v2 = [tx + sideLength, ty];

v3 = [tx + sideLength/2, ty + h_tri];

% Add Sierpinski triangle (recursive)

CSX = AddSierpinskiTriangle(CSX, 'Patch', v1, v2, v3, iter, h);

%% Ports (lumped)

start = [feed_x, feed_width/2, h];

stop  = [feed_x, feed_width/2, h + 0.5];

[CSX, port] = AddLumpedPort(CSX, 5, 1, [], start, stop, 50);

%% Mesh

mesh.x = SmoothMeshLines([0 tx + sideLength], max_res, 1.4);

mesh.y = SmoothMeshLines([0 ty + h_tri], max_res, 1.4);

mesh.z = SmoothMeshLines([0 h 10], max_res);

CSX = DefineRectGrid(CSX, unit, mesh);

%% Excitation and boundary

FDTD = InitFDTD();

FDTD = SetGaussExcite(FDTD, (f_start+f_stop)/2, f_stop-f_start);

FDTD = SetBoundaryCond(FDTD, {'PML_8','PML_8','PML_8','PML_8','PML_8','PML_8'});

%% Write and run

Sim_Path = 'SierpinskiMicrostrip';

Sim_CSX = 'sierpinski_microstrip.xml';

[~,~] = rmdir(Sim_Path, 's');

mkdir(Sim_Path);

WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);

RunOpenEMS(Sim_Path, Sim_CSX);

%% Post-process S11

port = calcPort(port, Sim_Path, f_start, f_stop, 1001);

figure;

plot(port.f/1e9, 20*log10(abs(port.uf.ref ./ port.uf.inc)));

xlabel('Frequency (GHz)'); ylabel('|S_{11}| (dB)');

title('S11 of Sierpinski Resonator with Microstrip Feed');

grid on;

%% --- Subfunction ---

function CSX = AddSierpinskiTriangle(CSX, metal_name, p1, p2, p3, iter, h)

    if iter == 0

        % AddPolygon expects 3D points

        CSX = AddPolygon(CSX, metal_name, 1, ...

            [p1, h; p2, h; p3, h]);

    else

        m12 = (p1 + p2)/2;

        m23 = (p2 + p3)/2;

        m31 = (p3 + p1)/2;

        CSX = AddSierpinskiTriangle(CSX, metal_name, p1, m12, m31, iter-1, h);

        CSX = AddSierpinskiTriangle(CSX, metal_name, m12, p2, m23, iter-1, h);

        CSX = AddSierpinskiTriangle(CSX, metal_name, m31, m23, p3, iter-1, h);

    end

end