function y = c_ext(x,epsilon,alpha,k,n,np)
ro=2*x*(n-1);
beta=atan(np/(n-1));
rp=4*x*np;

Qc=2-4*exp(-ro.*tan(beta)).*(cos(beta).*sin(ro-beta)./ro+(cos(beta)./ro).^2.*cos(ro-2*beta))+4*(cos(beta)./ro).^2.*cos(2*beta);
y = Qc.*pi.*x.^2./k.^2.*epsilon^4.*alpha.^5./(epsilon.^2-1)./x.^5.*sqrt((epsilon.^2-1)./(alpha^2*epsilon.^2-x.^2));
