function xk = waveno(wfreq,h)
% 
%    xk = waveno(wfreq,h);
%
% m-file function to calculate wavenumber
%
% inputs:  wfreq = array of wave frequencies (rad/sec) 
%	       h = mean water depth (m)
%
% output:     xk = array of resulting wave numbers (rad/m)
%
 
g = 9.81;   % acceleration of gravity (m/s^2)

% Calculate terms in wave speed approximation equation

y = h.*wfreq.*wfreq/g;
y1 = .66667 * y;
y2 = .35550 * y.^2;
y3 = .16084 * y.^3;
y4 = .06320 * y.^4;
y5 = .02174 * y.^5;
y6 = .00654 * y.^6;
y7 = .00171 * y.^7;
y8 = .00039 * y.^8;
y9 = .00011 * y.^9;

% Compute wave speed approximation, then wave number

c = sqrt(g*h./(y+1./(1+y1+y2+y3+y4+y5+y6+y7+y8+y9)));
xk = wfreq./c;


% Extra calculations (presently not included in function)
%xkh = xk .* h;
%xkh2 = xkh*2;