%
%  w07foil.m -- (djm:  13 feb 2002)
%    - with lift parameter & airfoil
%

%  initialize workspace

clear;  close all

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

%  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}
    circ = 1;
    Psi = zz + 1./zz + circ*(i*2*pi/2/pi)*log(zz);
    uu = (1 - 1./r2) + (2*yy.*yy)./(r2.^2) + circ*(2*yy./r2);
    vv =             - (2*xx.*yy)./(r2.^2) - circ*(2*xx./r2);

    ftitle = ['\bf flow around a cylinder'];
 
 case{2}
    lift = 0;
    cU = exp(-i*pi/6);  cUc = conj(cU);  cL = -2*pi*imag(cUc);
    Psi = (cU /2) *        (zz + sqrt(zz-1).*sqrt(zz+1)) ... 
	+ (cUc/2) *        (zz - sqrt(zz-1).*sqrt(zz+1)) ...
	- lift*(i*cL/2/pi) * log(zz + sqrt(zz-1).*sqrt(zz+1));

    ftitle = ['\bf flat airfoil'];
end

%  plot level curves of Psi

clf
contour(x,y,real(Psi),-8:0.4:5,'g--')
hold on
contour(x,y,imag(Psi),-5:.25:5)
colormap(autumn);  colorbar
contour(x,y,imag(Psi),[0 0],'k')

%  label plot

title(ftitle)
xlabel('\bf x-axis')
ylabel('\bf y-axis')
axis equal;
axis([-L L -L L])

if(flo==1 & circ~=0)
%    plot([-10 0],[0 0],'color','k','linewidth',4)
    text(-2.8,0.2,'\bf branch problem here!')
end

if(flo==2)
    plot([-1 1],[0 0],'color','k','linewidth',4)
    if(lift~=0)
        text(-2.9,0.2,'\bf branch problem here!')
        a = text(-0.5, 0.5,'\bf fast & lo');  set(a,'color','r')
        a = text(-0.5,-0.5,'\bf slow & hi');  set(a,'color','r')
    end
end