function alfa=SaderAlfafuncCil(rho,eta,omega,omega_0,rho_c,W,CN)
C1=1.8751041;
B=C1*sqrt(omega./omega_0).*sqrt(sqrt((1+((rho./rho_c)*GammaCilinder(rho,eta,omega,W))))); % I don't really know what B is physically, but it's necessary
alfa=ones(length(CN),length(omega)); %initialization
for jj=1:length(CN) %calculate alfa for one mode at a time
alfa(jj,:)= (2*sin(CN(jj)).*tan(CN(jj)))./(CN(jj).*((CN(jj))^4 -B.^4).*(sin(CN(jj))+sinh(CN(jj))));
end

