function [eta,vx,vy,ax,ay]=wvprop(theta,y,dol0,hohb,hol0,lol0,coef)
% ******************************************************************
% [eta,vx,vy,ax,ay]=wvprop(theta,y,dol0,hohb,hol0,lol0,coef)
% ************************************************************
%  Compute profile and wave kinematics on a "theta by y" grid
%     theta = 2*pi*((x/L)-(t/T)), (may be column vector or scalar)
%     y = vertical coordinate, measured vertically up from sea floor,
%         divided by deep water Airy wave length. (may be column
%         vector or scalar)
%     dol0 = depth over deep water Airy wave length
%     hol0 = wave height over deep water Airy wave length
%     lol0 = wave length over deep water Airy wave length
%     coef = stream function series coefficients, X(n)
%  Output:
%     eta = water level elevation above mean water level over wave height
%     vx = horizontal component of velocity times period over height
%     vy = vertical component of velocity times period over height
%     ax = horizontal component of acceleration times T^2 over H
%     ay = vertical component of acceleration times T^2 over H
%  ************************************************************

% Establish size of kinematics arrays
% ------------------------------------
[nr,nn]=size(y);
[nc,nn]=size(theta);

% Compute wave profile, eta/H
% ----------------------------
eta=profile(theta,dol0,hohb);

% Create y, theta, and mask matrices
% ----------------------------------
thet=pi*theta/180;
Y=zeros(nr,nc);
T=zeros(nr,nc);
mask=zeros(nr,nc);
for j=1:nc
   yeta=dol0+eta(j,1)*hol0;
   for i=1:nr
      if y(i,1)<=yeta
         mask(i,j)=1;
         Y(i,j)=2*pi*y(i,1)/lol0;
         T(i,j)=thet(j,1);
      end
   end
end

% Zero work arrays
% ------------------
vx=zeros(nr,nc);
vy=zeros(nr,nc);
CS=zeros(nr,nc);
SC=zeros(nr,nc);

% Compute series summations
% --------------------------
nt=numcoef(dol0,hohb)
for n=1:nt
   vx=vx+n*coef(n,1)*cosh(n*Y).*cos(n*T);
   vy=vy+n*coef(n,1)*sinh(n*Y).*sin(n*T);
   CS=CS+n*n*coef(n,1)*cosh(n*Y).*sin(n*T);
   SC=SC+n*n*coef(n,1)*sinh(n*Y).*cos(n*T);
end

% Calculate partial derivative matrices
% -------------------------------------
tpi=2*pi;
M1=tpi*tpi/lol0;
M2=M1*tpi/lol0;
pxx=M1*CS;
pyx=-M1*SC;
pxy=-M2*SC;
pyy=-M2*CS;

% Compute and mask kinematics arrays
% -----------------------------------
vx=-M1*vx;
vy=-M1*vy;
ax=tpi*((hol0/lol0)*vx-ones(nr,nc)).*pxx+hol0*vy.*pxy;
ay=tpi*((hol0/lol0)*vx-ones(nr,nc)).*pyx+hol0*vy.*pyy;
vx=vx.*mask;
vy=vy.*mask;
ax=ax.*mask;
ay=ay.*mask;