%
% bi-quad antenna
% 2.6GHz
%
% openEMS 64bit -- version v0.0.35-64-g3d684ba
% octave        -- version 5.2.0
% (C) 2020 Yellow Rabbit

close all
clear
clc

output_precision(5);

threads_num = 2;
%debug = ' --no-simulation';
%debug = ' --debug-PEC --no-simulation';
debug = ' ';
show = 1;

%% prepare simulation folder
Sim_Path = tilde_expand('~/tmp/b26');
Sim_CSX = 'bi_quad_ant.xml';

confirm_recursive_rmdir(0);
[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder

%% setup the simulation
physical_constants;
unit = 1e-3; % all length in mm

% frequency range of interest
f_start =  2.2e9;
f0      =  2.6e9;
f_stop  =  2*f0 - f_start;
freq = linspace(f_start, f_stop, 401); 
end_crit = 1e-6;
lambda = c0/f0/unit;

% wire
wh = 2.2;
ww = 1.3;

Rc = 8/2;
Ri = Rc - ww/2;
Ro = Rc + ww/2;
Rc_proj = Rc/sqrt(2);

port_length         = 1 + ww;
port_proj           = port_length * sqrt(2);
dR                  = Rc*sqrt(2) - Rc - port_length/2;
quad_size           = lambda/4 * 0.973;
delta_port          = (2*pi*Rc)/4;
quad_small          = (2*quad_size - dR*sqrt(2) - 2*delta_port)/2;
quad_big            = (2*quad_size + dR*sqrt(2) - 2*delta_port)/2;
reflector_off       = lambda/8.25 + wh/2;
reflector_thickness = wh; 
reflector_size      = quad_size*4 * 1.1;

quad_mesh           = ww/2;
port_mesh           = port_length/6;
max_res             = lambda/15;
reflector_mesh      = lambda/20;

Feed_R = 50;

% size of the simulation box
SimBox = [300 800 800];

%% setup FDTD parameter & excitation function
FDTD = InitFDTD('endCriteria', end_crit);
FDTD = SetGaussExcite(FDTD, 0.5*(f_start + f_stop), 0.5*(f_stop - f_start));
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; % boundary conditions
FDTD = SetBoundaryCond(FDTD, BC);

%% setup CSXCAD geometry & mesh
CSX = InitCSX();

% create fixed lines for the antenna port
mesh.x = [-wh -wh/2 0 wh/2 wh];
mesh.y = [-Rc_proj-ww/2 -Rc_proj+ww/2 0 Rc_proj-ww/2 Rc_proj+ww/2];
mesh.z = [-port_length/2-ww/2 -port_length/2+ww/2 0 port_length/2-ww/2 port_length/2+ww/2];
mesh = SmoothMesh(mesh, port_mesh, 1.3);

% create fixed lines for the antenna outline
mesh.x = [mesh.x -reflector_off-reflector_thickness -reflector_off];
mesh.y = [mesh.y -(3*Rc + quad_big + quad_small)/sqrt(2)+Rc_proj-Rc-ww/2 ...
                 -(3*Rc + quad_big + quad_small)/sqrt(2)+Rc_proj-Rc+ww/2 ...
                  (3*Rc + quad_big + quad_small)/sqrt(2)-Rc_proj+Rc-ww/2 ...
                  (3*Rc + quad_big + quad_small)/sqrt(2)-Rc_proj+Rc+ww/2];
mesh.z = [mesh.z -(Rc + quad_small)/sqrt(2)+Rc_proj-Rc-ww/2 ...
                 -(Rc + quad_small)/sqrt(2)+Rc_proj-Rc+ww/2 ...
                  (Rc + quad_small)/sqrt(2)-Rc_proj+Rc-ww/2 ...
                  (Rc + quad_small)/sqrt(2)-Rc_proj+Rc+ww/2];
mesh = SmoothMesh(mesh, quad_mesh, 1.3);

% create reflector mesh
mesh.y = [mesh.y -reflector_size/2 reflector_size/2];
mesh.z = [mesh.z -reflector_size/2 reflector_size/2];
mesh = SmoothMesh(mesh, reflector_mesh, 1.3);

% add air box
mesh.x = [mesh.x -SimBox(1)/4 SimBox(1)*3/4];
mesh.y = [mesh.y -SimBox(2)/2 SimBox(2)/2];
mesh.z = [mesh.z -SimBox(3)/2 SimBox(3)/2];

mesh = SmoothMesh(mesh, max_res, 1.5);

CSX = DefineRectGrid(CSX, unit, mesh);

% create a biquad
CSX = AddMaterial(CSX,'metal'); 
#CSX = AddMetal(CSX,'metal');
%% Copper
CSX = SetMaterialProperty( CSX, 'metal', 'Kappa', 56e6);
% inner
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 quad_big],
             'Transform', {'Rotate_X', -pi/4, 'Translate', ...
                 ['0,' num2str(Rc_proj) ',' num2str(Rc_proj - dR)]});
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 -quad_big],
             'Transform', {'Rotate_X',  pi/4, 'Translate', ...
                 ['0,' num2str(Rc_proj) ',' num2str(-Rc_proj + dR)]});
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 quad_big],
             'Transform', {'Rotate_X', pi/4, 'Translate', ...
                 ['0,' num2str(-Rc_proj) ',' num2str(Rc_proj -dR)]});
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 -quad_big],
             'Transform', {'Rotate_X',  -pi/4, 'Translate', ...
                 ['0,' num2str(-Rc_proj) ',' num2str(-Rc_proj + dR)]});
% outer
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 -quad_small],
             'Transform', {'Rotate_X',  pi/4, 'Translate', ...
              ['0,' num2str(3*Rc_proj + quad_big/sqrt(2)) ',' num2str((quad_small + Rc)/sqrt(2))]});
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 -quad_small],
             'Transform', {'Rotate_X', -pi/4, 'Translate', ...
              ['0,' num2str(-3*Rc_proj - quad_big/sqrt(2)) ',' num2str((quad_small + Rc)/sqrt(2))]});
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 quad_small],
             'Transform', {'Rotate_X', pi/4, 'Translate', ...
              ['0,' num2str(-3*Rc_proj - quad_big/sqrt(2)) ',' num2str(-(quad_small + Rc)/sqrt(2))]});
CSX = AddBox(CSX, 'metal', 10, [-wh/2 -ww/2 0], [wh/2 ww/2 quad_small],
             'Transform', {'Rotate_X',  -pi/4, 'Translate', ...
               ['0,' num2str(3*Rc_proj + quad_big/sqrt(2)) ',' num2str(-(quad_small + Rc)/sqrt(2))]});
% curve corners
p(1, 1) = Ri;
p(2, 1) = -wh/2;
p(1, 2) = Ri;
p(2, 2) = wh/2;
p(1, 3) = Ro;
p(2, 3) = wh/2;
p(1, 4) = Ro;
p(2, 4) = -wh/2;
CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [3*pi/4 5*pi/4], 'Transform', ...
               {'Translate', ['0,0,' num2str(Rc + port_length/2)]});
CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [7*pi/4 9*pi/4], 'Transform', ...
               {'Translate', ['0,0,' num2str(-Rc - port_length/2)]});
CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [5*pi/4 7*pi/4], 'Transform', ...
               {'Translate', ['0,' num2str((quad_big + quad_small + 2*Rc)/sqrt(2)) ',0']});
CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [pi/4 3*pi/4], 'Transform', ...
               {'Translate', ['0,' num2str(-(quad_big + quad_small + 2*Rc)/sqrt(2)) ',0']});

CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [7*pi/4 9*pi/4], 'Transform', ...
               {'Translate', ['0,' num2str(-(quad_big + 2*Rc)/sqrt(2)) ',' ...
                   num2str((Rc + quad_small)/sqrt(2) - Rc_proj)]});
CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [7*pi/4 9*pi/4], 'Transform', ...
               {'Translate', ['0,' num2str((quad_big + 2*Rc)/sqrt(2)) ',' ...
                   num2str((Rc + quad_small)/sqrt(2) - Rc_proj)]});
CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [3*pi/4 5*pi/4], 'Transform', ...
               {'Translate', ['0,' num2str(-(quad_big + 2*Rc)/sqrt(2)) ',' ...
                   num2str(-(Rc + quad_small)/sqrt(2) + Rc_proj)]});
CSX = AddRotPoly(CSX, 'metal', 10, 'y', p, 'x', [3*pi/4 5*pi/4], 'Transform', ...
               {'Translate', ['0,' num2str((quad_big + 2*Rc)/sqrt(2)) ',' ...
                   num2str(-(Rc + quad_small)/sqrt(2) + Rc_proj)]});
% reflector plane
CSX = AddMetal(CSX,'reflector'); 
start = [-reflector_off-reflector_thickness -reflector_size/2 -reflector_size/2];
stop = [-reflector_off reflector_size/2 reflector_size/2];
CSX = AddBox(CSX, 'reflector', 10, start, stop);

%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
start = [0 0 -port_length/2+ww/2];
stop  = [0 0  port_length/2-ww/2];
[CSX port] = AddLumpedPort(CSX, 10, 0, Feed_R, start, stop, [0 0 1], true);

%% nf2ff calc
start = [mesh.x(9) mesh.y(9) mesh.z(9)];
stop  = [mesh.x(end-8) mesh.y(end-8) mesh.z(end-8)];
[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);


%% write openEMS compatible xml-file
WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);


%% show the structure
if (show)
	CSXGeomPlot([Sim_Path '/' Sim_CSX], ['--export-polydata-vtk=' Sim_Path]);
end

%% run openEMS
Settings.LogFile = [Sim_Path '/openEMS.log'];
options = ['--engine=multithreaded --numThreads=' num2str(threads_num) debug];
RunOpenEMS(Sim_Path, Sim_CSX, options, Settings);

% Postprocessing & do the plots
freq = linspace(f_start, f_stop, 201);
port = calcPort(port, Sim_Path, freq);

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

%% swr & resonant
swr = (1 + abs(S11)) ./ (1 - abs(S11));
[swr_min idx] = min(swr);
[Xmin idx1] = min(abs(imag(Zin)));
swr_fr = freq(idx);
res_fr = freq(idx1);
Zr_swr = Zin(idx);
Zr_res = Zin(idx1);
disp('Minimum SWR frequency:');
disp(swr_fr);
disp('Resonant frequency (jX=0):');
disp(res_fr);
disp('Impedance at minimum SWR:');
disp(Zr_swr);
disp('Impedance at resonant frequency:');
disp(Zr_res);
disp('Minimum SWR:');
disp(swr_min);

% S11
plot(freq/1e6, 20*log10(abs(S11)));
grid on
title('Reflection coefficient');
xlabel('frequency f / MHz');
ylabel('S11 / dB');
print([Sim_Path '/s11.png']);
close;

% plot Zin
figure
plot(freq/1e6, real(Zin), 'k-', 'Linewidth', 2);
grid on
title('feed point impedance');
xlabel('frequency f / MHz');
ylabel('resistance Z_{in} / Ohm');
print([Sim_Path '/Zin-real.png']);
close;

figure
plot(freq/1e6, imag(Zin), 'r--', 'Linewidth', 2);
grid on
title('feed point impedance');
xlabel('frequency f / MHz');
ylabel('reactance Z_{in} / Ohm');
print([Sim_Path '/Zin-imag.png']);
close;

%% calculate 3D far field pattern
phiRange = -180:2.5:180;
thetaRange =  0:2.5:180;

nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, thetaRange*pi/180, phiRange*pi/180);

disp( ['directivity: Dmax = ' num2str(10*log10(nf2ff.Dmax)) ' dBi'] );

% plot far-field pattern with Matlab
%figure
%plotFF3D(nf2ff, 'logscale', -20)

%%
disp( 'Dumping far-field pattern to vtk (use Paraview to visualize)...' );
DumpFF2VTK([Sim_Path '/Bi_Quad_Pattern.vtk'], nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)) * nf2ff.Dmax, thetaRange, phiRange, 'scale', 0.05);

% vim: expandtab: sw=4 ts=4: