%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Flow Around a Joukowski Airfoil Using Complex Variables %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set the flow velocity at large distance if it is not yet set
if exist('U') ~= 1
  U=1
end

% set the angle of attack in degrees if it is not yet set
if exist('alp') ~= 1
  alp=0
end
% convert to radians
alpha=alp/180*pi;

% set radius of the cylinder if it is not yet set
if exist('r0') ~= 1
    r0=1
end

% set a value for the circulation if it is not yet set
if exist('Gamma') ~= 1
    Gamma=0
end

% location of the trailing edge on the cylinder in degrees
if exist('th_te') ~= 1
    th_te=0
end
% convert to radians
theta_te=th_te/180*pi;

% create mesh points around the cylinder in the z plane
[th,r]=meshgrid([-.001:5:360]*pi/180,[r0-.0001:.05:3]);
[x,y]=pol2cart(th,r);
z=x+i*y;

% if the circulation is nonzero, correct it to the right value
if Gamma ~= 0
    Gamma=4*pi*U*r0*sin(alpha+theta_te)
end

% compute the complex velocity potential at those points
F=U*(z*exp(-i*alpha)+r0^2*exp(i*alpha)./z)+i*Gamma/(2*pi)*log(z/r0);

% compute the real streamfunction
psi=imag(F);

% compute the conformally mapped physical coordinate zeta = xi + i eta
eps=r0*exp(-i*theta_te)-1;
zeta=(z-eps)+1./(z-eps);

% compute the corresponding real coordinates
xi=real(zeta);
eta=imag(zeta);

% plot the streamlines in the z plane
figure(1)
contour(x,y,psi,floor(min(min(psi))):.25:ceil(max(max(psi))))
axis("equal")
title("z-plane")

% plot the streamlines in the zeta-plane
figure(2)
contour(xi,eta,psi,floor(min(min(psi))):.25:ceil(max(max(psi))))
axis equal
title("zeta-plane")
