%
%  w07cmap.m -- (djm:  11 feb 2004)
%

%  initialize workspace
clear;  %close all

%  set x & y-vectors &  grid matrices ("help meshgrid")
L = 4.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;

%  flow around cylinder
Phi = (zz + 1./zz);

%  plot level curves of Phi in z-plane
figure(1);  clf
subplot(2,2,1);  hold on
axis([-1 1 -1 1]*L);  axis square
ylabel('\bf z-plane')
contour(x,y,real(Phi),[-2.0:0.5:2.0],'g--')
contour(x,y,imag(Phi),[-1.5:0.5:1.5],'k--')

th=0:2*pi/100:2*pi;
fill(cos(th),sin(th),'w')
plot(cos(th),sin(th),'m','linewidth',2)

plot(1,0,'rx')
plot(-1,0,'ro')
plot(0,1,'bx')
plot(0,-1,'bo')
plot(0,(1+sqrt(5))/2,'b*')

%  plot Phi-plane
subplot(2,2,2);  hold on
axis([-1 1 -1 1]*L);  axis square
ylabel('\bf \Phi-plane')
contour(x,y,real(zz),[-2.0:0.5:2.0],'g--')
contour(x,y,imag(zz),[-1.5:0.5:1.5],'k--')

contour(x,y,imag(zz),[0 0],'k')

%  plot boundary & points
plot([-2 2],[0 0],'m','linewidth',2)
plot(2,0,'rx')
plot(-2,0,'ro')
plot(0,0,'bx')
plot(0,0,'bo')
plot(0,1,'b*')

%  second example
%subplot(2,2,3);  hold on
figure(2);  clf;
subplot(2,1,1);  hold on
axis equal
axis([-3 3 0 3]);  
ylabel('\bf z-plane')

%  plot streamlines in z-plane
alpha = -11:0.2:5;
for psi = 0.25:0.5:5.25
    Phi = alpha + i*psi;
    z = (i*sqrt(1-Phi.^2) + log(Phi+i*sqrt(1-Phi.^2)) )/pi;
%    z = (i*sqrt(1-Phi.^2) + acosh(Phi) )/pi;
    plot(z,'b')
end

%  plot velocity potential in z-plane
beta = 0.001:0.2:6;
for phi = -10:0.5:4
    Phi = phi + i*beta;
%    z = (i*sqrt(1-Phi.^2) + log(Phi+i*sqrt(1-Phi.^2)) )/pi;
    z = (i*sqrt(1-Phi.^2) + acosh(Phi) )/pi;
    plot(z,'g--')
end

%  plot boundary & points
plot([-3 0 0 3],[1 1 0 0],'m','linewidth',2)
plot( 0,0,'rx')
plot( 0,1,'ro')
P0 = 4+i*0.0001;
%z0 = (i*sqrt(1-P0^2) + log(P0+i*sqrt(1-P0^2)) )/pi
z0 = (i*sqrt(1-P0^2) + acosh(P0) )/pi;
plot(z0,'bx')
P0 = -4+i*0.0001;
%z0 = (i*sqrt(1-P0^2) + log(P0+i*sqrt(1-P0^2)) )/pi
z0 = (i*sqrt(1-P0^2) + acosh(P0) )/pi;
plot(z0,'bo')

P0 = i*2.25;
z0 = (i*sqrt(1-P0^2) + acosh(P0) )/pi;
plot(z0,'k*')

%  plot contours in Phi-plane
subplot(2,1,2);  hold on

%  plot velocity potential pre-images in Phi-plane
beta = 0:0.2:6;
for phi = -10:0.5:4
    Phi = phi + i*beta;
    plot(Phi,'g--')
end

%  plot streamline pre-images in Phi-plane
alpha = -11:0.2:5;
for psi = 0.25:0.5:5.25
    Phi = alpha + i*psi;
    plot(Phi,'b')
end

axis equal
axis([-12 6 0 6])
ylabel('\bf \Phi-plane')

%  plot boundary & points
plot([-12 6],[0 0],'m','linewidth',2)
plot( 1,0,'rx')
plot(-1,0,'ro')
plot( 4,0,'bx')
plot(-4,0,'bo')
plot(0,2.25,'k*')