function [wk,err] = wavenumber(g,wa,wd,h,u,v,tol)
% Solves the wave dispersion relation sig^2 = g*k*tanh(k*h)
%
%Usage:
%  wk = wavenumber(g,wa,wd,h,u,v,tol)
%  [wk,err] = wavenumber(g,wa,wd,h,u,v)
%
%Description:
% Solves the wave dispersion relation 
%   sig^2 = g*wk*tanh(wk*h)
% where 
%   g = gravity [L/T^2]
%   h = water depth [L]
%   sig = Relative angular frequency [rad/T]
%   sig = wa - wk*cos(wd)*u - wk*cos(wd)*v = wa - wk*uk [rad/T]
%   uk = cos(wd)*u + sin(wd)*v [L/T]
%   u = current velocity in x direction [L/T]
%   v = current velocity in y direction [L/T]
%
% Input: 
%   wa = absolute angular frequency [rad/T]
%   wd = wave direction, going to [rad]
%   h = water depth [L]
%   u = current velocity in x direction [L/T]
%   v = current velocity in y direction [L/T]
%
% Output:
%   wk = wave number [rad/L]
%   err = error defined as sig^2 - g*k*tanh(k*h) [rad^2/T^2]
%
% The Newton-Raphson Method is given by
%   wk(n+1) = wk(n) - f(k(n))/fp(k(n))
% where
%   f = g*wk*tanh(wk*h) - sig^2
%
% Makes initial guess using 
% Guo, J. (2002) Simple and explicit solution of wave 
%   dispersion, Coastal Engineering, 46(2), 71-74.
%
% Author: Alex Sanchez, USACE-CHL
% Date: November 1, 2012
% Version 1.0

twopi = 2*pi;
%Input
if nargin==0
    g = 9.81; %gravity [m/s^2]
    T = 10; %wave period [s]
    freq = 1/T; %frequency [1/s]
    wa = twopi*freq; %absolute angular frequency [rad/s]
    wd = 0*pi/180; %wave angle [rad]
    h = 6; %water depth [m]
    u = -1;  %current velocity in x direction [m/s]
    v = 0;  %current velocity in y direction [m/s]
end
%Make initial guess with Guo (2002)
xi = wa/sqrt(g/h); %note: =h*wa/sqrt(h*g/h)
yi = xi*xi/(1.0-exp(-xi^2.4908))^0.4015;
wki = yi/h; %Initial guess without current-wave interaction
%current velocity in wave direction
uk = cos(wd)*u + sin(wd)*v; 
%Simple correction for currents (tends to overcorrect, hense the relaxation)
wki = 0.3*wki + 0.7*(wa-wki*uk)^2/g/tanh(wki*h);
%save for iteration
wk = wki; 
if(nargin==0)
    Li = twopi/wki %initial wave length [m]
end
if nargin<6
    tol = 0.0001; %Initalize loop variables
end
for iter=1:10
    sig = wa-wk*uk; 
    tanhkh = tanh(wk*h);
    f = g*wk*tanhkh - sig^2;
    fp = g*(tanhkh - h*wk*(tanhkh^2-1)) + 2*uk*sig;
    wk = wki - f/fp;
    if(wk<0 || wk>1e6) %No solution, assume wave blocking
        wk = -1; err = -1;
        break
    end
    wchg = abs(wk-wki);
    if(wchg<tol), break, end
    wki = wk;
end

if nargin==0
    L = twopi/wk   %final wave length [m]
    sig = wa-wk*uk;
    err = sig^2 - g*wk*tanh(wk*h) %Check
    iter
end

return

