%
%  code13.m -- (djm:  13 mar 2002)
%

%  initialize workspace & set constants

clear;  close all

L = pi;  ds = L/32;
H = L/16;  K = 1.0;  U = 1.0;

%  set coordinates & grid matrices

x = -L:ds:+L;  y = -4*H:ds:4*H;  t = 0;

[xg,yg] = meshgrid(x,y);

%  calculate amplitudes

B  = H*exp(+K*U*t); 
A1 = -(1+i)*U * B;
A0 =  (1-i)*U * B;

et =  B * exp(i*K*x); 

[eg,yg] = meshgrid(et,y);

%  calculate velocities (u = phi_x, v = phi_y)

u1g =  U + real(i*K *A1*exp(-K*yg).*exp(i*K*xg) ); 
v1g =      real(- K *A1*exp(-K*yg).*exp(i*K*xg) ); 

u0g = -U + real(i*K *A0*exp( K*yg).*exp(i*K*xg) ); 
v0g =      real(  K *A0*exp( K*yg).*exp(i*K*xg) ); 

%  plot velocity field

clf
quiver(x,y, u0g.*(eg>yg), v0g.*(eg>yg),0.5,'k')
hold on
quiver(x,y,-u0g.*(eg>yg),-v0g.*(eg>yg),0.5,'k.')
quiver(x,y, u1g.*(eg<yg), v1g.*(eg<yg),0.5,'k')
quiver(x,y,-u1g.*(eg<yg),-v1g.*(eg<yg),0.5,'k.')

axis equal;  axis([-L +L -4*H 4*H]);

plot(x,real(et),'r')

title('\bf kelvin-helmholtz instability (k=1, U=1)')
xlabel('\bf x-axis')
ylabel('\bf y-axis')

%v = v0g.*(eg<yg) + v1g.*(eg>yg);
%pcolor(x,y,v);  shading interp;  colorbar('horiz')