%  Translation flow around cylinder.
%	Domain is between theta = 0 and theta = pi for 1 < r < beta
clear all
clf reset
beta=input('beta = ')
figure(1)
clf reset
u=load('u.dat');
v=load('v.dat');
JMAX=size(u);
del=1/(JMAX(1)-1);
z=0:del:1;
gamma=1/log(beta);
[th,r]=meshgrid(0:del*pi:pi,exp(z./gamma));
[x,y]=pol2cart(th,r);
for i=1:JMAX
   for j=1:JMAX
      vr(i,j)=u(j,i);
      vth(i,j)=v(j,i);
      costh=cos(th(i,j));
      sinth=sin(th(i,j));
      vx(i,j)=vr(i,j)*costh-vth(i,j)*sinth;
      vy(i,j)=vr(i,j)*sinth+vth(i,j)*costh;
   end
end
quiver(x,y,vx,vy,'k')
hold on
plot(x(1,:)', y(1,:)','k')% Plot of boundary
axis image
xlabel('x')
ylabel('y')
title(['Flow Past Cylinder, beta = ',num2str(beta)])
pause 
%	View value of arrays
fname=input('Enter name of file to be viewed: ');
figure(2)
clf reset
array=load(fname);
contour(x,y,array',20,'k')
axis image
title(['Flow Past Cylinder, beta = ',num2str(beta),' Contours of ',num2str(fname)])
hold on
contour(x,y,array',[0 0],'r:')
plot(x(1,:)', y(1,:)','k')% Plot of boundary
