%
%  w06plate.m -- (djm:  14 feb 2004)
%    - useful stuff for hw #5
%

%  initialize workspace
clear;

%  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) problem A')
  disp('    2) problem B')
  disp(' ')

  flo = input(' choice:  ');
end

%  complex potential
switch(flo)

%  problem A
  case{1}
    Phi = -i*sqrt(i*zz-1).*sqrt(i*zz+1);
%    Phi = sqrt(zz.^2+1);
    ftitle = ['\bf problem A'];
%    cI = 40;  cR = 10;
    cI = -2:0.2:2;  cR = -2:0.2:2;
	
  case{2}
    Q = 2.0;
    Phi = zz + Q*(log(-zz)+i*pi)/(2*pi);
    ftitle = ['\bf problem B'];
%   cI = 40;  cR = 10;
    cI = -2:0.25:2;  cR = -5:0.5:5;
end

%  plot level curves of Phi
clf
contour(x,y,real(Phi),cR,'g--')
hold on
contour(x,y,imag(Phi),cI)

axis equal;  axis square
colormap(autumn);  colorbar

%  label plot
title(ftitle)
xlabel('\bf Re(z)')
ylabel('\bf Im(z)')

%  fancy stuff
if(flo==1)
    plot([0 0],[1 -1],'k','linewidth',2)
else
    contour(x,y,imag(Phi),[1 1]*(Q/2),'k')
end