%
%  code05.m -- (djm:  04 feb 2002)
%

%  initialize workspace

clear;  close all

L = 2.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;  rr2 = xx.^2 + yy.^2;

%  complex potential

flo = 6;

switch(flo)

%  flow in an eccentric annulus
  case{1}
    Phi = -(i/2/pi) * log((i*zz-1) ./ (i*zz+1));
    ftitle = ['\bf flow around an eccentric annulus'];
    co = 40;
	
%  flow through a gap
  case{2}
    Phi = log(zz + i*sqrt( -(zz.^2 - 1) ));
    ftitle = ['\bf flow through a gap'];
    co = 40;
	
%  flow around ellipse
  case{3}
    beta = pi/4;
    sq  = sqrt(zz-1).*sqrt(zz+1);
    Phi = exp(-i*beta) * (zz + sq) + exp( i*beta) * (zz - sq) *(1.5) ;
    ftitle = ['\bf flow around an ellipse'];
    co = 40;
	
%  flow over a post	
  case{4}
    Phi = i*sqrt(i*zz-1).*sqrt(i*zz+1);
    ftitle = ['\bf flow over a post'];
    co = 40;
	
  case{5}
    Phi = zz + 1./zz;
    ftitle = ['\bf flow past a circle'];
    co = 40;

  case{6}
    Phi = zz - 0.002*cos(pi*zz);
    ftitle = ['\bf wavy channel flow'];
    co = -1.5:0.1:1.5;

  case{7}
    Phi = zz.^(0.6);
    ftitle = ['\bf branch cut?'];
    co = 40;
end

%  plot level curves of Phi

clf
contour(x,y,real(Phi),10,'g--')
hold on

contour(x,y,imag(Phi),co,'b-')
%contour(x,y,imag(Phi),-1.5:.1:1.5,'b-')

contour(x,y,imag(Phi),[-1e-6 1e-6],'k-')
axis equal;  axis square

%  label plot

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