I = sqrt(-1);
gray = 210/255*ones(1,3);
theta = linspace(0,2*pi,100);

% %================================
% %Plot different discs in z-plane
% %================================
% figure('Color','w'); 
% line([-11 11], [0 0], 'Color','k'); hold on;
% line([0 0], [-11 11], 'Color', 'k');
% 
% rh = [1./(10:-2:2) 1 2:2:10];
% Th = pi/6*[1 2 4 5 7 8 10 11];
% 
% for i=1:length(rh)
%     z1 = rh(i)*cos(theta) + I*rh(i)*sin(theta);
%     if abs(rh(i) - 1) < 1e-12
%         plot(real(z1),imag(z1),'r', 'LineWidth',2);
%     else        
%         plot(real(z1),imag(z1),'b');
%     end
% end
% 
% for i=1:length(Th)
%     z1 = (1:0.2:10)*exp(I*Th(i));
%     plot(real(z1),imag(z1),'k');
% end
% 
% %plot(0,0,'ok', 'MarkerFaceColor','k');
% text(11.15,-0.1, 'Re(z)');
% text(-0.1, 11.15, 'Im(z)');
% axis([-11 11 -11 11]);
% axis square;
% axis off;
% 
% %==========================================
% %Plot the images of circles in the w-plane
% %==========================================
% 
% figure('Color','w'); 
% line([-6 6], [0 0], 'Color','k'); hold on;
% line([0 0], [-6 6], 'Color', 'k');
% 
% for i=1:length(rh)
%     z1 = rh(i)*cos(theta) + I*rh(i)*sin(theta);   
%     w1 = 0.5*( z1 + 1./z1);
%     
%     if abs(rh(i) - 1) < 1e-12
%         plot(real(w1),imag(w1),'r', 'LineWidth', 2); hold on;
%     else        
%         plot(real(w1),imag(w1),'b');
%     end
% end
% 
% for i=1:length(Th)
%     z1 = (1:0.2:10)*exp(I*Th(i));
%     w1 = 0.5*( z1 + 1./z1);   
%     
%     plot(real(w1),imag(w1),'k');
% end
% 
% %plot(0,0,'ok', 'MarkerFaceColor','k');
% text(6.2,-0.2, 'Re(w)');
% text(0.2, 6.2, 'Im(w)');
% axis([-6 6 -6 6]); 
% axis off;


% ==========================================================================
% ==========================================================================

%==============================
%Plot the Upper Half--plane
%==============================
%x = [-10 10 10 -10];
%y = [-10 -10 10 10];
figure('Color','w'); 
%fill(x,y,gray, 'EdgeColor',gray); hold on;
line([-10 10], [0 0], 'Color', 'k'); hold on;
line([0 0],[-10 10], 'Color', 'k');
plot([-1,1],[0,0],'r','LineWidth',2);

xh = -10:0.05:10;
ch = [-logspace(0.5,-1.5,21) logspace(-1.5,0.5,21)];
for i=1:length(ch)
    yh1 = ch(i)*ones(size(xh));
    plot(xh,yh1,'b');
end

text(4.2,-0.2, 'Re(z)');
text(-0.2, 2.2, 'Im(z)');
axis([-4 4 -2 2]); 
axis off;

%===================================
% Plot the Upper half plane, 
% with the cylinder 
%===================================
z0 = -0.2*(1-I);
r0 = 0.2*sqrt(37);

figure('Color', 'w'); 
t = linspace(0,2*pi,50);
x = cos(t); y = sin(t);
zb = x + I*y;
zb = z0 + r0*zb;
zb = 0.5*( zb + 1./zb);

x_t = real(zb); y_t = imag(zb);

fill(x_t, y_t, gray, 'EdgeColor', 'r'); hold on;
line([-10 10], [0 0], 'Color', 'k'); 
line([0 0],[-10 10], 'Color', 'k');

xh_p = logspace(-6,1,200);
xh_n = -logspace(1,-6,200);
ch = [-logspace(0.5,-1.5,21) logspace(-1.5,0.5,21)];
for i=1:length(ch)
    yh1 = ch(i);
    w1 = xh_p+I*yh1;
    w2 = xh_n+I*yh1;
    
    %============================================
    % Map straight lines to lines around cylinder
    %============================================
    z1 = w1+sqrt(w1.^2-1);
    z2 = w2-sqrt(w2.^2-1);
    
    %==================================================
    % Map 0-centerd cylinder to off centered cylinder 
    %==================================================
    z1 = z0 + r0*z1;
    z2 = z0 + r0*z2;    

    %=====================================================
    % Map the off-centered cylinder to Joukowski airfoil
    %=====================================================
    z1 = 0.5*(z1 + 1./z1);
    z2 = 0.5*(z2 + 1./z2);
        
    plot(real(z1),imag(z1),'b');
    plot(real(z2),imag(z2),'b');
end

text(2.1,-0.2, 'Re(z)');
text(-0.2, 2.1, 'Im(z)');
axis([-2 2 -2 2]); 
axis off;

