%
%  w08drift.m -- (djm:  23 feb 2004)
%

%  initialize workspace & set constants

clear;  close all

global H K sp

L = 2*pi;  ds = L/20;  N = 13;
H = 0.2;  K = -1.0;  G = 1.0;
sp = sqrt(G/abs(K));

%  initialize code parameters & initial particles

x = -L/2:ds:+L/2;
T0 = 0.0;  dt = 2*pi/N;

Lpt = [-L/2:L/4:L/2];
Zi = [Lpt + i*H*cos(K*(Lpt+sp*0*dt)) ; Lpt-i/2 ; Lpt-i];

%  initial plot
eta = H*cos(K*(x-sp*0*dt));
wave = plot(x,eta,'b--');  hold on;  
plot(real(Zi),imag(Zi),'r*');
pts = plot(real(Zi),imag(Zi),'r*');
axis equal;  axis([-L/2 +L/2 -L/4 3*H])

%  label plot
title(['\bf particle trajectories & stokes drift!!   (h=0.2, k=1.0, speed=1.0)'])
xlabel('\bf x-axis')
ylabel('\bf y-axis')

disp(' ');  disp('  return to continue')
disp(' ');  pause;  disp('  thanks')
figure(1)

%  time loop (ODE for particle paths)
for j=1:3*N
  [T,Z] = ode45('w08drift_ode', [T0 T0+dt], Zi);
  oldZi = Zi;
  Zi = Z(end,:);  T0 = T0 + dt;
  
%  update plot
  set(wave,'erase','xor');
  eta = H*cos(K*(x+sp*j*dt));
  set(wave,'ydata',eta);

  set(pts,'erase','xor');
  plot(oldZi,'k.','markersize',7)
  set(pts,'xdata',real(Zi),'ydata',imag(Zi));  
  pause(1)
end

plot(Zi,'g*')