%  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

%  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 = 1.;
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:10:150;
nx = max(size(x));
z = 0:-10:-70;
nz = max(size(z));

dt = 0.3;

for jt = 1:90
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

figure(1)
clf reset
quiver(xq, zq, u, w)
axis([0 150 -70 10])
ylabel('depth, m')
hold on
plot(x, 5.*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(10, -58, dep)
timtex = ['time, sec = ', num2str(time)];
text(10, 20, timtex)
title('Fluid velocity and surface displacement of a gravity wave')
axis('equal')

M(jt) = getframe;

end

movie(M, 2, 6)