[kx,ky]=kaes(f,h,z,Nmod,Nit,delta,rho,c)

function [kx,ky]=kaes(f,h,z,Nmod,Nit,delta,rho,c)

% Esta funcion obtiene los Nmod primeros modos kx y ky para el caso del canal
% en dos dimensiones con altura h e impedancia en la parte superior z, la parte
% inferior se supone rigida. la frecuencia es f, Nit es el numero de iteraciones 
% maximo al buscar la solucion de los modos usando el metodo de Newton y delta
% es el valor de variacion de la solucion con la cual se detiene la iteracion.


if (abs(z)<0.001)
    L=1:(2*Nmod);
    W=L*pi-pi/2;
else
    beta=2*pi*f*h*rho/z;
    L=1;
    if (imag(beta)>=0)
        W(L)=beta;
        dw=1;
        N=0;
        while ((N<Nit) & (abs(dw)>delta))
            auxnum=j*(W(L)*(exp(j*2*W(L))-1)+beta*(exp(j*2*W(L))+1));
            auxden=(W(L)*(exp(j*2*W(L))+1)+(beta-j)*(exp(j*2*W(L))-1));
            dw=auxnum/auxden;
            W(L)=W(L)+dw;
            N=N+1;
        end;

       W(L)=abs(real(W(L)))+abs(imag(W(L)))*j;
       aux=W(L)*sinkom(W(L))-j*beta*coskom(W(L));
        if (aux<0.00001)
            L=L+1;
        end;
    end;
    M=0;
    while (L<=(2*Nmod))
        if (abs(beta)<5000 & (M>0) & (L<=2*Nmod))
            W(L)=M*pi+j*beta/(M*pi);
            dw=1;
            N=0;
            while ((N<Nit) & (abs(dw)>delta))
                auxnum=j*(W(L)*(exp(j*2*W(L))-1)+beta*(exp(j*2*W(L))+1));
                auxden=(W(L)*(exp(j*2*W(L))+1)+(beta-j)*(exp(j*2*W(L))-1));
                dw=auxnum/auxden;
                W(L)=W(L)+dw;
                N=N+1;
            end;

            W(L)=abs(real(W(L)))+abs(imag(W(L)))*j;
            aux=W(L)*sinkom(W(L))-j*beta*coskom(W(L));
            if (aux<0.00001) & (length(find(abs(W(L)-W)<0.1))<=1)
                L=L+1;
            end;
        end;
        if (beta~=j & (L<=2*Nmod))
            W(L)=M*pi+pi/2-(pi/2+M*pi)/(1+j*beta);
            dw=1;
            N=0;
            while ((N<Nit) & (abs(dw)>delta))
                auxnum=j*(W(L)*(exp(j*2*W(L))-1)+beta*(exp(j*2*W(L))+1));
                auxden=(W(L)*(exp(j*2*W(L))+1)+(beta-j)*(exp(j*2*W(L))-1));
                dw=auxnum/auxden;
                W(L)=W(L)+dw;
                N=N+1;
            end;
%           
                 W(L)=abs(real(W(L)))+abs(imag(W(L)))*j;
              
            aux=W(L)*sinkom(W(L))-j*beta*coskom(W(L));
            if (aux<0.00001) & (length(find(abs(W(L)-W)<0.1))<=1)
                L=L+1;
            end;
        end;
        M=M+1;
    end;
end;
kyaux=W/h;
kxaux=sqrt((2*pi*f/c)^2-(kyaux).^2);
kxaux(find(imag(kxaux)>0))=-kxaux(find(imag(kxaux)>0));
[kaux,IX]=sort(abs(imag(kxaux)));
kxaux=kxaux(IX);
kyaux=kyaux(IX);
kx=kxaux(1:Nmod).';
ky=kyaux(1:Nmod).';
k=[kx,ky];