function k=kfromw(w,p)
%KFROMW  Wavenumber as a function of omega for surface waves
%   K = KFROMW(W,P) gives the wavenumber K (in m^-1) for surface waves
%   of pulsation W (in s^-1), using the parameter structure P. By default,
%   P is set for an air-water interface with infinite depth (see
%   WAVE_PARAMETER for the default values). If K is a vector, W is a vector
%   of same length.
%
%   Example:
%     f = linspace(10,50);
%     lambda = 2*pi./kfromw(f*2*pi);
%     plot(f,lambda*1000);
%     xlabel('f (Hz)');
%     ylabel('\lambda (mm)');
%     axis([0 50 0 25])
%
%   See also WAVE_PARAMETER, OMEGA, KCKG
%
%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 1.21,  Date: 2010/04/28

% 2009/02/10: v1.00
% 2010/04/15: v1.10 : works also for finite depth (p.h0)
% 2010/04/20: v1.20 : use wave_parameter
% 2010/04/28: v1.21 : pre-allocation of k



error(nargchk(1,2,nargin));

if nargin==1
    p=[];
end

p = wave_parameter(p);   % fills the default values in p

k = zeros(1,numel(w)); % pre-allocate

if isinf(p.h0)  % Case of infinite depth: solve the cubic polynom  
    for i=1:numel(w)
        % solves the cubic equation
        ro = roots([p.gamma/p.rho, 0, p.g, -w(i)^2]);
        k(i) = ro(3);  % take the third root (real root)
    end
    
else    % Case of finite depth: find the zero of w-w(k)
    % first, find the solution for infinite depth
    kinf = kfromw(w,wave_parameter(p,'h0',inf));
    
    % second, find the zero of w-w(k), starting from the infinite
    % depth solution
    for i=1:numel(w)
        k(i) = fzero(@(kk) w(i) - sqrt(tanh(kk*p.h0) * (p.g*abs(kk) + p.gamma*abs(kk)^3/p.rho)), kinf(i));
    end
end
