ap=1
as=30
Os=2.35
epsilon=sqrt(10^(ap/10)-1)
k1=sqrt(10^(ap/10)-1)/sqrt(10^(as/10)-1)
n=ceil(acosh(1/k1)/acosh(Os))
k=1/cosh(acosh(1/k1)/n)
B=(sqrt(1+(epsilon/k1)^2)+(epsilon/k1))^(1/n)
a=(B-1/B)/2
b=(B+1/B)/2
%format long
mu=(1:2:2*n-1)'
p=1./(-a*sin(mu*pi/(2*n))+i*b*cos(mu*pi/(2*n)))/k
Hjme=real(poly(p))             %jmenovatel
if rem(n,2)                    %n - liche
    mu=mu([1:end]~=(n+1)/2)
    z=i./cos(mu*pi/(2*n))/k    %nuly pren. funkce
    Hcit=real(poly(z))*n*k1/epsilon/k %citatel pren. funkce
    Pcit=eye(1,n)*k/n/k1       %citatel char. funkce
else                           %n - sude
    z=i./cos(mu*pi/(2*n))/k    %nuly pren. funkce
    Hcit=real(poly(z))/sqrt(epsilon^2/k1^2+1) %citatel pren. funkce
    Pcit=eye(1,n)/k1           %citatel char. funkce
end
Pjme=real(poly(z))             %jmenovatel char. funkce