function eta = eta(p,x,t)
% eta Return wave amplitude as a function of position and time.
m = length(x(:,1));


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);


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


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

% wavep = 2*0.001*(1/(8*pi))*1025*(9.81^2)*trapz(f,S./f);
% disp(['Available Wave Power = ',num2str(wavep),' kW'])

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

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



eta = 0.5*H(1)*cos(k(1)*x - p.beta*2*pi*f(1)*t + p.phases(1));


for i = 2:length(f)

eta = eta +  0.5*H(i)*cos(k(i)*x - p.beta*2*pi*f(i)*t + p.phases(i));

end

%for j = 1:length(x)
%eta(j) = eta(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
eta = p.A*((1-exp(-t./tau)).*cos(p.beta*w*t)*cos(k*x)+(1-exp(-t./tau)).*sin(p.beta*w*t)*sin(k*x));
%eta = p.A*(cos(w*t)*cos(k*x)+ sin(w*t)*sin(k*x));

case 2  %Create timeseries based upon supplied 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; %wavenumber from Dispersion relation


%for j = 1:length(t);
     eta = zs(1)*cos(k(1)*x - p.beta*2*pi*p.freq(1)*t + p.phases(1));
  
     for i = 2:length(p.freq)
          %disp(['i= ',num2str(i),'  j= ',num2str(j)]);
          eta = eta +  zs(i)*cos(k(i)*x - p.beta*2*pi*p.freq(i)*t + p.phases(i));
     end
     
%end    


otherwise
error('Invalid IncWave type')       
        
end