function [w,c,cg,p]=omega(k,p)
%OMEGA  Dispersion relation for surface waves
%   W = OMEGA(K,P) gives the pulsation W (in s^-1) for surface wave
%   of wavenumber K (in m^-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.
%
%   [W,C,CG,P] = OMEGA(K,P) also returns the phase velocity C=W/K, the group
%   velocity CG=dW/dK, and the parameter structure P.
%
%   If gravity (p.g) is negative, complex pulsation is returned. The
%   imaginary part is the growth rate of the instabililty.
%
%   Example with default parameters:
%     k = logspace(2,3,200);
%     c = omega(k)./k;
%     plot(2*pi./k, c);
%     xlabel('Wavelength (m)');
%     ylabel('Phase velocity (m/s)');
%
%   Example with user-specified parameters:
%     k = logspace(2,3,200);
%     p = wave_parameter('gamma',0.06,'h0',0.1);
%     c = omega(k,p)./k;
%     plot(2*pi./k, c);
%
%   See also WAVE_PARAMETER, KFROMW, KCKG

%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 1.21,  Date: 2009/02/28

% 2007/10/18: v1.00
% 2008/07/02: v1.10, added phase and group velocities
% 2009/02/10: v1.11, small bug in help text
% 2010/04/20: v1.20, use wave_parameter
% 2010/04/28: v1.21, cosmetics


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

if nargin==1
    p=[];
end

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

% dispersion relation for infinite depth:
w = sqrt(p.g * abs(k) + p.gamma * abs(k).^3 / p.rho);

% finite depth correction
if ~isinf(p.h0)
    w = w .* sqrt(tanh(k*p.h0));
end

% additional output parameters: phase and group velocities
if nargout>1
    c = w./k; % phase velocity
    eps = 1e-4; % small parameter for finite difference
    w2 = omega((1+eps)*k,p);
    w1 = omega((1-eps)*k,p);
    cg = (w2-w1)./(2*eps*k); % group velocity
end
