ap=1
as=30
Os=2.35
%format long
epsilon=sqrt(10^(ap/10)-1)
k1=sqrt(10^(ap/10)-1)/sqrt(10^(as/10)-1)
k=1/Os
n=ceil(ellipke(1-k1^2)/ellipke(k1^2)*ellipke(k^2)/ellipke(1-k^2))
a=ellipke(1-k1^2)/ellipke(k1^2)
options = optimset('TolX',1e-10);
k=fminbnd(@(x) (ellipke(x^2)*a-ellipke(1-x^2)*n)^2,0,1,options)
Kk=ellipke(k^2)
if rem(n,2)
    m=(n-1)/2
    mu=1:m
    Ozm=ellipj(2*mu/n*Kk,k^2)
    k1=k^n*prod((1-Ozm.^2)./(1-k^2*Ozm.^2))^2
else
    m=n/2
    mu=1:m
    Ozm=ellipj((2*mu-1)/n*Kk,k^2)
    k1=k^n*prod(Ozm)^4
end
Opm=1./(k*Ozm)
asn=10*log10(1+epsilon^2/k1^2)
vinf=ellipke(k^2)/ellipke(k1^2)/n *ellipticf(atan(1/epsilon),1-k1^2)
[SN,CN,DN]=ellipj(vinf,1-k^2);
alpha0=-SN/CN
alpha=alpha0*sqrt(1-Ozm.^2).*sqrt(1-k^2*Ozm.^2)./(1+alpha0^2*k^2*Ozm.^2)
beta=Ozm*sqrt(1+alpha0^2)*sqrt(1+k^2*alpha0^2)./(1+alpha0^2*k^2*Ozm.^2)
if rem(n,2)
    H0=-alpha0*prod(alpha.^2+beta.^2)/prod(Opm)^2
    Hjme=poly([alpha+i*beta,alpha-i*beta,alpha0])/H0
    Pcit=Hjme(1)*poly([i*Ozm,-i*Ozm,0])
else
    H0=1/sqrt(1+epsilon^2)*prod(alpha.^2+beta.^2)/prod(Opm)^2
    Hjme=poly([alpha+i*beta,alpha-i*beta])/H0
    Pcit=sqrt(Hjme(1)^2-1)*poly([i*Ozm,-i*Ozm])
end
Hcit=poly([i*Opm,-i*Opm])
Pjme=Hcit