% (1) Display Joukowski transform of a circle.

% Unit circle
thetas = 0:(pi/300):2*pi;
zeta_x = cos(thetas);
zeta_y = sin(thetas);
zeta = zeta_x + i.*zeta_y;

% Joukowski transform of unit circle is a flat plate.
z = zeta + 1./zeta;
z_x = real(z);
z_y = imag(z);

% Uncambered airfoil
zeta2_x = 1.1*zeta_x - 0.08;
zeta2_y = 1.1*zeta_y;
zeta2 = zeta2_x + i.*zeta2_y;
z2 = zeta2 + 1./zeta2;
z2_x = real(z2);
z2_y = imag(z2);

% Cambered airfoil
M = 1.12; % magnification of the cylinder for the cambered airfoil.
delta_z = -0.1 + 0.15*i; % offset of the cylinder for the cambered airfoil.
zeta3_x = M*zeta_x + real(delta_z);
zeta3_y = M*zeta_y + imag(delta_z);
zeta3 = zeta3_x + i.*zeta3_y;
z3 = zeta3 + 1./zeta3;
z3_x = real(z3);
z3_y = imag(z3);


figure(1)

hold off
plot(zeta_x, zeta_y, 'k:', 'LineWidth', 2)
hold on
plot(z_x, z_y, 'k', 'LineWidth', 2)

plot(zeta2_x, zeta2_y, 'b:', 'LineWidth', 2)
plot(z2_x, z2_y, 'b', 'LineWidth', 2)

plot(zeta3_x, zeta3_y, 'r:', 'LineWidth', 2)
plot(z3_x, z3_y, 'r', 'LineWidth', 2)

lgnd=legend('Unit Circle','Flat Plate','1.1x-0.08, 1.1y','Uncambered Airfoil','           1.12x-0.1, 1.12y+0.15','Cambered Airfoil');
legend('boxoff')
set(lgnd,'FontSize',12)
set(lgnd,'Position',[0.39,0.45,0.25,0.1])

set(gca,"FontSize",16,"LineWidth",1,"TickDir",'in')
axis equal

print('joukowski1.pdf','-dpdf','-tight','-S200,150','-F:Helvetica:2.5')



% (2) Calculate and display the 2D flow field around a cylinder. This is
% the offset, oversized cylinder that will be transformed later into the airfoil.

aoa = pi/18; % angle of attack
v0 = 1.0;    % background speed
Gamma = 3.8; % circulation

x1 =  -4:0.02:4;
y1 = (-2:0.02:2)';
Nx = numel(x1);
Ny = numel(y1);
x_arr = ones(Ny,1) * x1;
y_arr = y1 * ones(1,Nx);
z_arr = x_arr + i*y_arr;
% Recall the cylinder offset used for the cambered foil in figure 1:
%  zeta3_x = M*zeta_x + real(delta_z);
%  zeta3_y = M*zeta_y + imag(delta_z);
% We now invert that coordinate change as follows...
zC_arr = (z_arr - delta_z)./M;
% ... and use the new zC coordinate to do the flow around the cylinder:
w_cylinder = v0.*(zC_arr*exp(-i*aoa)+1./(zC_arr*exp(-i*aoa))) + i*Gamma*log(zC_arr)/(2*pi) ;

figure(2)
hold off
contour(x1, y1, real(w_cylinder), -6:0.1:6, 'g', 'LineWidth', 1)
hold on
contour(x1, y1, imag(w_cylinder), -2:0.1:4, 'b', 'LineWidth', 1)
plot(zeta3_x, zeta3_y, 'r', 'LineWidth', 2) % overplot cylinder from figure 1
axis equal
% set(gca,"FontSize",12,"LineWidth",2,"TickDir",'out',"TickLength",[0.015,0.015])
set(gca,"FontSize",16,"LineWidth",1,"TickDir",'in')
print('joukowski2.pdf','-dpdf','-tight','-S250,150','-F:Helvetica:2.5')


% (3) Implement the Joukowski transform, z = zeta + 1/zeta, yielding the flow around 
% an airfoil. Unfortunately, we need the inverse transform, which isn't unique. (I'm
% forced into this because I need to plot the complex potential as a function of x,y
% rather than zeta_x, zeta_y.) The key was in realizing which root of the quadratic
% is needed where, and making an empirical formula to do it automatically.

which_root = sign(x_arr); % Almost right
pesky_wedge = (x_arr > -0.7) & (x_arr < 2) & (y_arr >= 0) & (y_arr < 0.1*(x_arr+0.7).*(2-x_arr));
which_root = (1 - 2*pesky_wedge) .* which_root; %improved
zeta_arr = ( z_arr + which_root .* sqrt( z_arr.^2 - 4 ) )./2;
figure(4)
contour(x1,y1,real(zeta_arr),-5:0.1:5)
figure(5)
contour(x1,y1,imag(zeta_arr),-5:0.1:5)
% Recall the cylinder offset used for the cambered foil in figure 1:
%  zeta3_x = M*zeta_x + real(delta_z);
%  zeta3_y = M*zeta_y + imag(delta_z);
% We now invert that coordinate change as follows...
zetaC_arr = (zeta_arr - delta_z)./M;
% So that...
w_foil = v0.*(zetaC_arr*exp(-i*aoa)+1./(zetaC_arr*exp(-i*aoa))) + i*Gamma*log(zetaC_arr)/(2*pi) ;

figure(3)
hold off
contour(x1, y1, real(w_foil), -6:0.1:6, 'g', 'LineWidth', 1)
hold on
contour(x1, y1, imag(w_foil), -2:0.1:4, 'b', 'LineWidth', 1)
plot(z3_x, z3_y, 'r', 'LineWidth', 2) % overplot airfoil outline for comparison
axis equal
% set(gca,"FontSize",12,"LineWidth",2,"TickDir",'out',"TickLength",[0.015,0.015])
set(gca,"FontSize",16,"LineWidth",1,"TickDir",'in')
print('joukowski3.pdf','-dpdf','-tight','-S250,150','-F:Helvetica:2.5')




