function [k, p, u, v] = eqwave(m, omega, y, r_or_l)
% function [k, p, u, v] = eqwave(m, omega, y, r_or_l)
%
% Calculate the wavenumber and the meridional variation of
% pressure and velocity for an equatorial wave, given the
% meridional mode number and the frequency.  All quantities
% are nondimensionalized.
%
% Inputs:
%  m        meridional mode number; -1 for K wave, 0 for Y wave, etc.
%  omega    frequency
%  y        vector of latitudes
%  r_or_l   'r' for the RH side of the dispersion curve, 'l' for the left
%              (optional; default is 'r')
%
% Outputs:
%  k        wavenumber
%  p, u, v  pressure, zonal velocity, meridional velocity;
%              vectors corresponding to latitudes y.
%              p and u are real, v is imaginary.
%
%  If only the wavenumber is desired, then only m and omega are
%  required as inputs.
%
%  Nondimensionalization assumes length scale sqrt(c/beta),
%  time scale 1/sqrt(c*beta).
%

if nargin < 2 | (nargin == 2 & nargout > 1)
   help eqwave
   return
end
if nargin < 4,
   r_or_l = 'r';
end
if nargin < 3,
   y = 0;
end

[kp, km] = eqwave_k(m, omega);
if r_or_l(1) == 'r', k = kp; else k = km; end
if m == -1,
   v = 0;
   u = exp(-0.5 * y.^2);
   p = u;
elseif m == 0,
   u = hermfunc(1, y);
   p = u;
   v = (sqrt(2) / (i * omega)) * hermfunc(0, y);
else
   a1 = sqrt(2*(m+1)) * hermfunc(m+1, y) / (omega - k);
   a2 = sqrt(2*m) * hermfunc(m-1, y) / (omega + k);
   u = a1 + a2;
   p = a1 - a2;
   v = (2/i) * hermfunc(m, y);
end
