%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Plot the potential flow solution for the streamlines %%%%%%%%
%%%%%%%%%%% along an impulsively started circular cylinder %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This program assumes that an array of r-values r-vals has been set
% previously by another program.  (Normally, cylpsi1.m or cylpsi2.m)

% velocity at infinite distance
U=1;

% cylinder radius
r_0=1;

% create theta-values
th_vals=[-.001:5:360]*pi/180;

% create 2D mesh points
[th,r]=meshgrid(th_vals,r_vals);
[x,y]=pol2cart(th,r);

% flow velocity at the wall
u_e=2*U*sin(th_vals);

% compute the streamfunction profile
y_vals=r_vals-1;
psi0=.5*(r_vals-r_0^2./r_vals);

% very slightly modify the wall value to capture the wall streamline
psi0(1)=-.001*psi0(2);

% compute the streamfunction (note that this is a column-row vector
% multiplication producing a two-dimensional grid of psi-values.)
psi=psi0*u_e;

% plot the streamlines, with more densely spaced streamlines near zero
contour(x,y,psi,...
[floor(min(min(psi))):.1:-.1 -.005 0 .005 .1:.1:ceil(max(max(psi)))])
