function [deta_dx] = deta_dx(p,x,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(t);
if (max(dim) > 1)
    error('ERROR:  t must be a scaler')
end


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

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

deta_dx =  -H(1)*k(1)*sin(k(1)*x - beta*2*pi*p.freq(1)*t + phases(1));

for i = 2:length(p.freq)
deta_dx =  deta_dx-H(i)*k(i)*sin(k(i)*x - beta*2*pi*p.freq(i)*t + phases(i));
end

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