function [velx, vely, velz] = vel(p,x,z,t)
% vel Return velocity due to incident wave as a function of position and time.
dim = size(x);
if (min(dim) > 1)
    error('ERROR:  x must be a scaler or array')
end
if (dim(1) > 1)
    x = x';
end

dim = size(z);
if (min(dim) > 1)
    error('ERROR:  z must be a scaler or array')
end
if (dim(2) > 1)
    z = z';
end


dim = size(t);
if (max(dim) > 1)
    error('ERROR:  t must be a scaler')
end

H = (2*p.spectrum*p.df).^0.5;
k = (2*pi*p.freq).^2/p.g;

phases = p.phases;
beta = p.beta;

velx = H(1)*(2*pi*p.freq(1))*exp(k(1)*z)*cos(k(1)*x - beta*2*pi*p.freq(1)*t + phases(1));
vely = zeros(length(z),length(x));
velz = H(1)*(2*pi*p.freq(1))*exp(k(1)*z)*sin(k(1)*x - beta*2*pi*p.freq(1)*t + phases(1));

for i = 2:length(f)
  velx = velx + H(i)*(2*pi*p.freq(i))*exp(k(i)*z)*cos(k(i)*x - beta*2*pi*p.freq(i)*t + phases(i));
  velz = velz + H(i)*(2*pi*p.freq(i))*exp(k(i)*z)*sin(k(i)*x - beta*2*pi*p.freq(i)*t + phases(i));
end

if(p.tau > 0)
  for j = 1:length(x)
    velx(j) = velx(j).*(1-exp(-t(j)/p.tau));
    velz(j) = velz(j).*(1-exp(-t(j)/p.tau));
  end
end
              
end