%
%  w04flow.m - velocity field (djm: 28 jan 04)
%
%

%  initialize workspace & set constants (try different values of B)
clear;  close all

L = pi;  ds = 2*L/45;
kk = 1.0;  ll = 1.0;  K = sqrt(kk^2 + ll^2);

%  set x & y-vectors &  grid matrices ("help meshgrid")
x = -L:ds:+L;  y = -L:ds:+L;
[xx,yy] = meshgrid(x,y);
rr2 = xx.^2 + yy.^2;

%  calculate psi
psi = cos(kk*xx).*cos(ll*yy);

%  calculate velocities (u = psi_y, v = -psi_x)
uu = -ll*cos(kk*xx).*sin(ll*yy);
vv =  kk*sin(kk*xx).*cos(ll*yy);

%  plot vector field ("help quiver") outside unit circle ("help gt")
clf
quiver(x,y, uu, vv,0.3)
hold on
quiver(x,y,-uu,-vv,0.3,'.')
axis([-L +L -L +L]);  axis square;

%  plot level curves of vorticity (two versions, try both)
contour(x,y,K*psi,[-1:0.2:1]);  colorbar

%  plot one special contour
contour(x,y,K*psi,[0 0],'k');  colormap(jet)

%  label plot
title(['\bf 2D, steady, rotational flow'])
xlabel('\bf x-axis')
ylabel('\bf y-axis')