function K = gravd(W,h)
%GRAVD Gravity wave wavenumber determination
% k = rad/meter
% k = GRAVD(w,h)
% w= ang freq (can be a vector), h = water depth (meters)

if (length(W) < 1)
	error('You must specify at least one frequency');
end
G = 9.79329;
N = length(W);
W2 = W.*W;  kDEEP = W2/G;  kSHAL = W/sqrt(G*h);
erDEEP = ones(size(W)) - G*kDEEP.*dtanh(kDEEP*h)./W2;
one = ones(size(W));
d = one; done = zeros(size(W));	nd = find(done == 0);
k1 = kDEEP;  k2 = kSHAL; e1 = erDEEP; ktemp = done; e2 = done;
while 1 > 0
    e2(nd) =  one(nd) - G * k2(nd) .* dtanh(k2(nd)*h)./W2(nd);
    d = e2.*e2;
	done = d < 1e-20;
    if done
        K = k2; break;
    end
	nd = find(done == 0);
    ktemp(nd) = k1(nd) - e1(nd).*(k2(nd)-k1(nd))./(e2(nd)-e1(nd));
    k1(nd)=k2(nd);  k2(nd) = ktemp(nd);  e1(nd) = e2(nd);       % Reset values
end

function y = dtanh(x)
%DTANH	DTANH(X) is the more stable hyperbolic tangent of x.

a = exp(x.*(x<=50));	% Make sure nothing blows up to NaN
one = ones(size(x));
y = (abs(x) > 50) .* (abs(x)./x) + (abs(x) <= 50).*((a-one./a)./(a+one./a));
