%
%  w01flow.m - velocity field (djm: 07 jan 04)
%
%  save file as "w01flow.m" (make sure that it is on the matlab path!)
%  to run, type "w01flow" at the matlab prompt ">>"
%

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

L = 4.0;  ds = 2*L/45;
B = 0.0

%  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 = yy .* (1 - 1./(rr2)) + (B/2)*log(rr2);

%  calculate velocities (u = psi_y, v = -psi_x)
uu = (1 - 1./(rr2)) + (2*yy.*yy ./ (rr2.^2) ) + (B*yy ./ rr2 );
vv =                - (2*xx.*yy ./ (rr2.^2) ) - (B*xx ./ rr2 );

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

%  plot unit circle
th = 0:pi/50:2*pi;
plot(cos(th),sin(th),'k')

%  plot level curves of psi (two versions, try both)
contour(x,y,psi .* (rr2>=1),[-8.5:1:8.5]);  colorbar
%contour(x,y,psi,[-5.5:1:5]);  colorbar

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

%  label plot
title(['\bf my first vector plot (B=' num2str(B) ')'])
xlabel('\bf x-axis')
ylabel('\bf y-axis')

%  plot a point
xs = 0;  ys = 0;
plot(xs,ys,'r*')