%  Translational potential flow around cylinder.
%	Domain is between theta = 0 and theta = pi for 1 < r < beta
clear all
clf reset
beta=2;
u=load('u.dat');
v=load('v.dat');
JMAX=size(u);
% JMAX=[21 21];  %overide
del=1/(JMAX(1)-1);
[th,r]=meshgrid(0:del*pi:pi,1:(beta-1)*del:beta);
[x,y]=pol2cart(th,r);
vx=x;    % initialize vx and vy arrays
vy=x;
s=x;
vx=(1-cos(2.*th).*r.^(-2));
vy=-sin(2.*th).*r.^(-2);
quiver(x,y,vx,vy,'k')
hold on
plot(x(1,:)', y(1,:)','k')% Plot of boundary
axis image
xlabel('x')
ylabel('y')
title(['Potential Flow Past Cylinder, beta = ',num2str(beta)])
pause
% Plot stream function
s=y.*(1-1./(x.^2+y.^2));
clf reset
contour(x,y,s,20,'k')
title(['Potential Flow Past Cylinder, beta = ',num2str(beta), ...
       ' Contours of Stream Function'])
hold on
plot(x(1,:)', y(1,:)','k')% Plot of boundary
axis image

