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


switch p.Type
    
case 0  %Pierson-Moskowitz Spectrum
f = p.df:p.df:p.fmax;

alpha = .0081;
beta = 0.74;

f0= 1/p.T;% peak frequency


Hs = 2*p.A;
f0 = ((0.161*p.g/Hs)^0.5)/(2*pi);


%peak frequency from Significant Wave Height 
 
S = (alpha*p.g^2/(2*pi)^4)*exp(-beta*(f0./f).^4)./f.^5;

H = 2*(2*S*p.df).^0.5;


k = (2*pi*f).^2/p.g;


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

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

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

tau  = 2*p.T;  %Ramp up amplitude over two periods

for j = 1:length(x)

    velx(j) = velx(j).*(1-exp(-t(j)/tau));
    velz(j) = velz(j).*(1-exp(-t(j)/tau));

end

                
case 1  %Ramp to Monochromatic
w = 2*pi/p.T;
k = w^2/p.g;
tau  = 2*p.T;  %Ramp up amplitude over two periods

beta = p.beta;
velx = (1-exp(-t./tau))*exp(k*z)*p.A*w*(cos(k*x)*cos(beta*w*t)+sin(k*x)*sin(beta*w*t));
vely = zeros(length(z),length(x));
velz = (1-exp(-t./tau))*exp(k*z)*p.A*w*(sin(k*x)*cos(beta*w*t)-cos(k*x)*sin(beta*w*t));

case 2  %Use specified spectrum

df = p.freq(2)-p.freq(1);
zs = (2*p.spectrum*df).^0.5; %spectral energy to eta (wave amplitude)
k = (2*pi*p.freq).^2/p.g;

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

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


otherwise
error('Invalid IncWave type')       
        
end