% Calcuation of dimensions for a waveguide multi-hole directional coupler
% Based on 'The Design of Topwall Waveguide Directional Couplers',
% K. M. Hervey, DSTO Australia, research note ERL-0559-RN
%
% This implementation (C) 2018,2019 Tobias Plüss <tpluess@ieee.org>
% with thanks to Bernd Burchardt


clear all; close all; clc

% coupling in db
cc = 10;

% lower frequency
fl = 26.5e9;

% upper frequency
fh = 40e9;

% waveguide dimensions
a = 7.11e-3;
b = 3.55e-3;

% wall thickness
t = 0.2e-3;

% number of holes
N = 12;

% number of arrays
NA = 6;

% number of rows
NR = 2;

% distance of the holes to the walls, as fraction of a
k = 0.203;

% wavelengths at the two corner frequencies
c = 2.998e8;

% mean frequency (shouldn't this be the geometric mean?) and wavelength
f = (fl + fh) / 2;

% free space wavelength
lambda = @(f) c/f;

% guided wavelength
lg = @(f) 1/sqrt(1/lambda(f)^2 - 1/(2*a)^2);

x = k*a;

% lower and upper guided wavelength
ll = lg(fl);
lh = lg(fh);




ca = 1-c;


% calculate hole spacing.
sp = ll*lh/(2*(ll+lh));
fprintf('hole spacing: %f mm\n', sp*1000);

% electrical distance of the holes
phi = 2*pi*sp/ll;
ww = cos(phi);

aet = @(d) 1.0103*t + 0.0579*d/2;
amt = @(d) 1.0064*t + 0.0819*d/2;

ke = @(d,f) exp(-2*pi*aet(d)./(1.305*d) .* sqrt(1-(1.305*d/lambda(f)).^2));
km = @(d,f) exp(-2*pi*amt(d)./(1.705*d) .* sqrt(1-(1.705*d/lambda(f)).^2));

qe = @(d,f) 2*lambda(f)./(pi*1.305*d) .* tan(pi*1.305*d./(2*lambda(f)));
qm = @(d,f) 2*lambda(f)./(pi*1.705*d) .* tan(pi*1.705*d./(2*lambda(f)));

u1 = @(d,f) sin(pi*x/a)^2/6 .* km(d,f) .* qm(d,f);
u2 = @(d,f) cos(pi*x/a)^2*lg(f)^2/(24*a^2) .* km(d,f) .* qm(d,f);
u3 = @(d,f) -sin(pi*x/a)^2*lg(f)^2/(12*lambda(f)^2) .* ke(d,f) .* qe(d,f);
U = @(d,f) 2*pi*d.^3/(a*b*lg(f)) .* (u1(d,f) + u2(d,f) + u3(d,f));



% find the coupler polynomial.
poly = couplerpoly(N);
A = zeros(size(poly,1), 1);
for row = 1:size(poly,1)
    A(row) = polyval(poly(row,:), ww);
end



Aa = 20*log10(A / sum(A));

N2 = ceil(N/2);

% total number of holes per row
NH = N + (NA-1)*N2;
fprintf('number of holes per row: %d\n', NH);

AS = zeros(NA, NH);

for j=1:NA
    for k=1+(j-1)*N2 : NH
        ind = k-(j-1)*N2;
        if ind > length(A)
            AS(j,k) = 0;
        else
            AS(j,k) = A(k-(j-1)*N2);
        end
    end
end

s = sum(sum(AS));

if NA > 1
    co = 20*log10(s ./ sum(AS)) + cc + 6.0206*(NR-1);
else
    co = 20*log10(sum(AS) ./ AS) + cc + 6.0206*(NR-1);
end

d_initial = nthroot(12*a^2*b./(pi*10.^(co/20)),3);

func = @(d) norm(20*log10(1./U(d,f)) - co);
d_best = fminsearch(func, d_initial, optimset('MaxFunEvals', 30000, 'MaxIter', 30000));

% best and rounded hole diam.
d = round(50e3*d_best)/50e3;
fprintf('initial, best and rounded hole diameters:\n');
disp([d_initial' d_best' d'])

cnew = 20*log10(1./U(d,f));

% actual coupling
coup = 20*log10(sum(10.^(-cnew/20)))+6.0206*(NR-1);
fprintf('actual coupling: %f\n', coup);
fprintf('distance of the rows: %f mm\n', (a-2*x)*1e3)
fprintf('total length of coupled section: %f\n', NH*sp*1e3)

% simulation
f = linspace(fl, fh);

k = zeros(length(f), 4);


dist = (0:NH-1)*sp;

for ind = 1:length(f)



phase1 = -2*pi/lg(f(ind))*dist;
isol_ideal = 20*log10(abs(sum(U(d_best, f(ind)) .* exp(-1i*2*phase1))));
isol = 20*log10(abs(sum(U(d, f(ind)) .* exp(-1i*2*phase1))));

isol_max = 20*log10(abs(sum(U(d+0.02e-3, f(ind)) .* exp(-1i*2*phase1))));
isol_min = 20*log10(abs(sum(U(d-0.02e-3, f(ind)) .* exp(-1i*2*phase1))));

phase2 = fliplr(phase1);
phase = phase1 + phase2;
num = (0:NH-1);
koppel_ideal = 20*log10(abs(sum(NR*(U(d_best, f(ind)) .* exp(-2i *phi)))));
koppel = 20*log10(abs(sum(NR*(U(d, f(ind)) .* exp(-2i *phi)))));

koppel_max = 20*log10(abs(sum(NR*(U(d+0.02e-3, f(ind)) .* exp(-2i *phi)))));
koppel_min = 20*log10(abs(sum(NR*(U(d-0.02e-3, f(ind)) .* exp(-2i *phi)))));



k(ind, :) = [isol_ideal isol koppel_ideal koppel];


end

figure(1);
plot(f/1e9, k);
grid on
legend('isolation ideal', 'isolation rounded', 'coupling ideal', 'coupling rounded')
xlabel('f / GHz');
ylabel('dB');



montecarlo = 0;

xplus = 0.05e-3;
xminus = -0.05e-3;

r = [0 1; 1 1] \ [xplus; xminus];

res_isol = zeros(montecarlo, length(f));
res_koppel = res_isol;

for run = 1:montecarlo

    tol = r(1)*rand(size(d)) + r(2);



    for ind = 1:length(f)

        phase1 = -2*pi/lg(f(ind))*dist;
        isol = 20*log10(abs(sum(U(d+tol, f(ind)) .* exp(-1i*2*phase1))));

        phase2 = fliplr(phase1);
        phase = phase1 + phase2;
        num = (0:NH-1);

        koppel = 20*log10(abs(sum(NR*(U(d+tol, f(ind)) .* exp(-2i *phi)))));



        res_isol(run, ind) = isol;
        res_koppel(run, ind) = koppel;

    end


figure(2);

plot(f, res_isol, 'r', f, res_koppel, 'b');

end




if 0
op = fopen('output.txt', 'w');

for ind = 1:montecarlo


    for j = 1:length(f)

        fprintf(op, '%f %f %f\n', f(j), res_isol(ind, j), res_koppel(ind, j));
    end
    fprintf(op, '\n\n');


end

fclose(op);

end




