%
%  w06poten.m -- (djm:  03 feb 2004)
%

%  initialize workspace
clear;  %close all

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

%  menu
flo = [];  
while (isempty(flo)==1)
  disp(' ')
  disp('  flow choices:')
  disp('    1) around an annulus')
  disp('    2) through a gap')
  disp('    3) past an ellipse')
  disp('    4) over a post')
  disp('    5) past a circle')
  disp('    6) in a wavy channel')
  disp('    7) with a branch cut')
  disp('    8) source')
  disp('    9) image sources')
  disp(' ')

  flo = input(' choice:  ');
end

%  complex potential
switch(flo)

%  flow in an eccentric annulus
  case{1}
    Phi = -(i/2/pi) * log((i*zz-1) ./ (i*zz+1));
    ftitle = ['flow around an eccentric annulus'];
    cI = 40;  cR = 10;
	
%  flow through a gap
  case{2}
    Phi = log(zz + i*sqrt( -(zz.^2 - 1) ));
    ftitle = ['flow through a gap'];
    cI = 40;  cR = 10;
	
%  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 = ['flow around an ellipse'];
    cI = 40;  cR = 10;
	
%  flow over a post	
  case{4}
    Phi = i*sqrt(i*zz-1).*sqrt(i*zz+1);
    ftitle = ['flow over a post'];
    cI = 40;  cR = 10;
	
  case{5}
    Phi = zz + 1./zz;
    ftitle = ['\bf flow past a circle'];
    cI = -2:0.25:2;  cR = -5:0.5:5;

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

  case{7}
    Phi = zz + log(zz);
    ftitle = ['\bf branch cut?'];
    cI = 40;  cR = 10;
    
  case{8}
    Phi = log(zz)/(2*pi);
    ftitle = ['\bf source/sink'];
    cI = 40;  cR = 10;
    
  case{9}
    Phi = log((zz-1).*(zz+1)) / (2*pi);
    ftitle = ['\bf image sources'];
    cI = 40;  cR = 10;
end

%  plot level curves of Phi

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

contour(x,y,imag(Phi),cI)
contour(x,y,imag(Phi),[0 0],'k-')

axis equal;  axis square
colormap(autumn);  colorbar

%  label plot

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