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

%  initialize workspace & set constants

clear;  close all

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

%  set coordinates & grid matrices

x = -L:ds:+L;  y = -L/2:ds:H+ds;
[xg,yg] = meshgrid(x,y);

%  calculate velocities (u = psi_y, v = -psi_x)

ug = -K*H*sp* exp(K*yg).* cos(K*xg); 
vg = -K*H*sp* exp(K*yg).* sin(K*xg);

%  plot velocity field

clf
subplot(2,1,1)
quiver(x,y, ug, vg,0.5)
hold on
quiver(x,y,-ug,-vg,0.5,'.')

axis equal;  axis([-L +L -L/2 +3*H]);

%  plot wave surface & level curves of psi
%    denser coordinates & grid matrices

x1 = -L:ds/16:+L;  y1 = -L/2:ds/16:H+ds;
[xg1,yg1] = meshgrid(x1,y1);

eta =  H*                 cos(K*x1 ); 
psi =  H*sp* exp(K*yg1).* cos(K*xg1); 

plot(x1,eta,'k')
contour(x1,y1,psi,10)

%  label plot

title(['\bf left-moving water wave (h=0.2, k=1, speed=-1)'])
%xlabel('\bf x-axis')
ylabel('\bf y-axis')

%  masked plot

subplot(2,1,2)
quiver(x,y, ug, vg,0.5)
hold on
quiver(x,y,-ug,-vg,0.5,'.')

axis equal;  axis([-L +L -L/2 +3*H]);

%  plot wave surface & level curves of psi

sky = [-L x1 +L -L; 3*H eta 3*H 3*H];

plot(x1,eta,'k')
contour(x1,y1,psi,10)
fill(sky(1,:),sky(2,:),'w')

eta =  H*cos(K*(x1+pi/6)); 
plot(x1,eta,'k--')

%  label plot

title(['\bf left-moving wave with translation (dashed)'])
xlabel('\bf x-axis')
ylabel('\bf y-axis')