%
%  code06.m -- (djm:  11 feb 2002)
%

%  initialize workspace

clear;  close all

L = 3.0;  ds = 2*L/101;

%  set x & y-vectors &  grid matrices ("help meshgrid")

x = -L:ds:+L;  y = -L:ds:+L;
[xx,yy] = meshgrid(x,y);  
zz = xx + i*yy;
[tt rr] = cart2pol(xx,yy);  r2 = rr.^2;

%  complex potential

flo = 1;

switch(flo)

%  flow around a cylinder (circ = -2*pi)
  case{1}
    Phi = zz + 1./zz + (i*2*pi/2/pi)*log(zz);
    uu = (1 - 1./r2) + (2*yy.*yy)./(r2.^2) + (2*yy./r2);
    vv =             - (2*xx.*yy)./(r2.^2) - (2*xx./r2);

    ftitle = ['\bf flow around a cylinder (vel pot = green, strmfn' ...
	      ' = blue)'];
end

%  plot level curves of Phi

clf
contour(x,y,real(Phi),-2*pi:pi/8:pi,'g--')
hold on
contour(x,y,imag(Phi),-2:.2:4,'b-')
contour(x,y,imag(Phi),[0 0],'k-')

%  label plot

title(ftitle)
xlabel('\bf x-axis')
ylabel('\bf y-axis')

Nj = 4;  NN = Nj*floor(size(x,2)/Nj);
[xp,yp] = meshgrid(x(1:Nj:NN),y(1:Nj:NN));
rp2 = xp.^2 + yp.^2;
up = uu(1:Nj:NN,1:Nj:NN).*(rp2>1);
vp = vv(1:Nj:NN,1:Nj:NN).*(rp2>1);  

quiver(xp,yp,up,vp,0.3,'r');  quiver(xp,yp,-up,-vp,0.3,'r.');
axis equal;
axis([-L L -L L])

theta = 0:pi/50:2*pi;
fill(cos(theta),sin(theta),'w')

text(-2.5,0.2,'\bf branch problem here!')
text(-1.5,-1.5,'\bf \psi=0 separating streamline & 2 stagnation pts')