function p = couplerpoly(N)

% calculate the polynomials for a Chebyshev directional coupler
% the algorithm is described in https://hb9fsx.ch/files/chebyshev-coupler.pdf
% (C) 2018,2019 Tobias Plüss <tpluess@ieee.org>


if mod(N, 2)
    A = zeros((N+1)/2, (N+1)/2);
else
    A = zeros(N/2, N/2);
end

% go through all the relevant chebyshev polynomials;
% for example, for N=9 this would be T8,T6,T4,T2,T0
for ind = 1:N/2
    num = N-2*ind+1;
    col = [zeros(1, N-num-1) 2*chebypoly(num)]';
    A(:, ind) = col(1:2:end);
end
if mod(N, 2)
    row = [zeros(1, N-1) chebypoly(0)];
    A(end, :) =  A(end, :) + row(1:2:end);
end

p = chebypoly(N-1);
B = diag(p(1:2:end));

result = A \ B;

polys = zeros(size(result, 1), N);
polys(:, 1:2:end) = result(:, :);

numerator = polys(1, :);

halfmat = zeros(size(polys));

for ind = 1:size(polys, 1)
    denominator = [polys(ind, :) zeros(1, size(polys, 2)-1)];

    [q, ~] = deconv(denominator, numerator);

    halfmat(ind, :) = q;
end

if mod(N,2)
    halfrows = floor(N/2);
    result = [halfmat; flipud(halfmat(1:halfrows, :))];
else
    result = [halfmat; flipud(halfmat)];
end

p = fliplr(result);

end

%[q,r] = deconv(result(end,:), result(1,:))


% else
%     % N is even
% end
% 2*chebypoly(4)
% 2*[0 0 chebypoly(2)]
% 2*[0 0 0 0 chebypoly(0)]
