%	Irrotational, solenoidal flow in 2-D in the complex plane
%	w(z)=-(U-iV)z-(U+iV)a^2/z, Translation flow around cylinder.
%	Domain is between theta = 0 and theta = 2pi for r => a
%  Plot domain is a < r < b
clear all
clf reset
a=1;
b=3;
U=1;
V=1;
[th,r]=meshgrid((0:10:360)*pi/180,a:.02:b);
[x,y]=pol2cart(th,r);
z=x+i*y;
w=z; %allocating storage
theta=angle(z);
tsize=size(theta);
for j=1:(tsize(1)*tsize(2))
   if theta(j)<0 
      theta(j)=2*pi+theta(j);
   end
end
w=-(U-i*V).*z -a^2*(U+i*V).*(real(z)-i.*imag(z))./(z.*conj(z));
contour(x,y,real(w),29,'r')
hold on
contour(x,y,imag(w),15,'k')
plot(x(1,:)', y(1,:)','k')% Plot of boundary
axis image
xlabel('x')
ylabel('y')
title('Potentials and Streamlines for Flow Past Cylinder')
pause
p=-(-U+((a./r).^2).*(U.*cos(-2.*th)-V.*sin(-2.*th))).^2 ...
  -(+V+((a./r).^2).*(U.*sin(-2.*th)+V.*cos(-2.*th))).^2;
pcolor(x,y,p)
colormap(hot)
shading flat
axis image
