close all
clear all

A = .5;
T = 12;

%Compute incident wave profile
wave = IncWave;
wave = set(wave,'A',A);         %Amplitude = 1 meter
wave = set(wave,'T',T);
wave = set(wave,'df',.001);
wave = set(wave,'Type',3);      %Type 3 = Modified Pierson Moskowitz
wave = set(wave,'beta',-1);     %Waves moving towards negative x
wave = set(wave,'tau',0);     %Ramp waves up over two periods.

f = get(wave,'freq');
H = get(wave,'spectrum');
figure
plot(f,H)
xlabel('freq (Hz)');
ylabel('m^2/s');
title(['PM:  Hs =' num2str(2*A) '  Tp = ' num2str(T) ' s']);


dt = .1;
disp('Computing incident wave profile');
t = 0:dt:100;
eta0 = eta(wave,0,t);

figure
plot(t,eta0)
xlabel('time (s)')
ylabel('Wave Height at Origin (m)');
title(['PM:  Hs =' num2str(2*A) '  Tp = ' num2str(T) ' s']);


epsilon = 0.005;

x = -100:.1:100;
figure('Position',[   129         659        1435         289]);
pl1 = plot(x,NaN(size(x)));
ylim([-.5 .5]);
grid on;
xlabel('m');
ylabel('m');
title(['PM:  Hs =' num2str(2*A) '  Tp = ' num2str(T) ' s' '   t = 0s']);

hold all
pl2 = plot(1,NaN,'x');
for i = 1:200
eta0 = eta(wave,x,t(i));

deta_dx = diff(eta0)./diff(x);
d2eta_dx2 = diff(deta_dx)./diff(x(1:end-1));

idx = find((abs(deta_dx(1:end-1)) < epsilon) & (d2eta_dx2 < 0) & (eta0(1:end-2) > 0));


set(pl2,'XData',x(idx));
set(pl2,'YData',eta0(idx));

set(pl1,'YData',eta0);

title(['PM:  Hs =' num2str(2*A) '  Tp = ' num2str(T) ' s' '   t = ' num2str(i*dt) 's']);

pause(.1);

end