%  gwavemovie;  make a movie of wave fluid velocity 
%  A pmovie of a progating gravity wave.  wavelength =  
%  100 m;  you must input a water depth.   The surface
%  displacement is shown as a red line;  its amplitude
%  is not shown to the scale of the water depth or of the
%  particle velocity (i.e, amplitude is arbitrary).
%
%  Jim Price, January 2001
%

clear
clear memory
close all
 path(path, 'G:\educ\rotation')
 
%  Set some default graphics things.
set(0,'DefaultTextFontSize',14)
set(0,'DefaultAxesFontSize',14)
set(0,'DefaultAxesLineWidth',1.6)
set(0,'DefaultLineLineWidth',1.4)


g = 9.8;
denw = 1025.;
sig = 7.3e-2;                  %  surface tension, N/m
a = 3.;
lambda = 100. ;       %  wavelength, m
k = 2*pi./lambda;

H = input('  Enter the water depth (10 - 5000 m) ');

omega = sqrt((g*k + (sig/denw)*(k^3))*tanh(k*H));
T = 2*pi/omega;
c = sqrt((g/k + (sig/denw)*k)*tanh(k*H));

x = 0:5:100;
nx = max(size(x));
z = 0:-10:-70;
nz = max(size(z));


nt = 125;
xp = zeros(4,nt+1);
zp = xp;

xp(:,1) = [40. 40. 40. 40.];
zp(:,1) = [-5. -15. -25. -35.];


dt = 0.3;

for jt = 1:nt
time = jt*dt;
pht = omega*time;

eta = a*cos(x*k - pht);

for jx=1:nx
for jz=1:nz
u(jz,jx) = a*omega*(cosh(k*(z(jz) + H))/sinh(k*H))*cos(k*x(jx)-pht);
w(jz,jx) = a*omega*(sinh(k*(z(jz) + H))/sinh(k*H))*sin(k*x(jx)-pht);
if -z(jz) > H
u(jz,jx) = 0.; w(jz, jx) = 0.;
end
xq(jz,jx) = x(jx);
zq(jz,jx) = z(jz);
end
end

%  track a particle
for m=1:4
xph(m) = xp(m,jt) + 0.5*dt*a*omega*(cosh(k*(zp(m,jt) + H))/sinh(k*H))*cos(k*xp(m,jt)-pht);
zph(m) = zp(m,jt) + 0.5*dt*a*omega*(sinh(k*(zp(m,jt) + H))/sinh(k*H))*sin(k*xp(m,jt)-pht);

xp(m,jt+1) = xp(m,jt) + dt*a*omega*(cosh(k*(zph(m) + H))/sinh(k*H))*cos(k*xph(m)-(pht+dt/2));
zp(m,jt+1) = zp(m,jt) + dt*a*omega*(sinh(k*(zph(m) + H))/sinh(k*H))*sin(k*xph(m)-(pht+dt/2));
end

figure(1)
clf reset
set(gcf,'position', [100 100 600 500])
quiver(xq, zq, u, w, 0)
axis([0 80 -60 20])
ylabel('depth, m')
xlabel('x distance, m')
hold on
plot(x, eta, 'r')
hold on; plot([0 150], [0 0], ':')
hstr = num2str(H);
kh = floor(k*H*10)/10;
khstr = num2str(kh);
dep = ['\lambda = 100 m,   h = ', hstr, ' m,   kh = ', khstr ];
text(15, -55, dep)
timtex = ['time = ', num2str(time), ' sec'];
text(10, 10, timtex)
title('Eulerian velocity and particle trajectories')

for m=1:4
hold on
plot([xp(m,jt+1)], [zp(m,jt+1)], 'g*')
hold on
plot(xp(m,1:jt), zp(m,1:jt), 'r')
end

M(jt) = getframe(gcf, [0 0 600 500]);

end

%  print -depsc2 gwaveposter.eps

%  mpgwrite(M, colorbar, 'gwavemovie.mpg')