%
%  code10.m
%

%  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/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));
plot(x,eta,'b--');  hold on;  plot(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')

%  time loop (ODE for particle paths)

for j=1:3*N
  [T,Z] = ode45('code10a', [T0 T0+dt], Zi);
  Zi = Z(end,:);  T0 = T0 + dt;
  
%  update plot

  plot(x,eta,'w--');
  eta = H*cos(K*(x+sp*j*dt));
  plot(x,eta,'b--');  plot(Zi,'r*')

  pause(1)
  plot(Zi,'w*');  plot(Zi,'k.','markersize',7)
end

plot(Zi,'g*')