%
%  code02.m
%
%  save file as "code02.m" (make sure that it is on the matlab path!)
%  to run, type "code02" at the matlab prompt ">>"
%

%  initialize workspace & set constants(try B=0,1,2,3)

clear;  close all

L = 4.0;  ds = 2*L/45;
B = 1.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 (huh?)

contour(x,y,psi .* (rr2>=1),10)

%  label plot

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

%  plot a point

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