%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Plot the two-term Blasius expansion for the streamlines %%%%%%
%%%%%%%%%%% along an impulsively started circular cylinder %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% velocity at infinite distance
U=1;

% cylinder radius
r_0=1;

% kinematic viscosity
nu=0.001;

% time
t=1;

% initial conditions
tau_0=2/sqrt(pi);
tau_1=2/sqrt(pi)*(1+4/(3*pi));

% maximum eta value
eta_max=10;

% integrate the ODE
[eta_vals, zeta_vals] = ode45(@zetasys, [0 eta_max], [0;0;tau_0;0;0;tau_1]);

% create r- and theta-values
r_vals=r_0+sqrt(4*nu*t)*eta_vals;
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 just above the boundary layer and its x-derivative
% note the sign change from the fact that th_vals=180-theta
u_e=2*U*sin(th_vals);
u_ep=-2*U*cos(th_vals)/r_0;

% compute the streamfunction profiles
y_vals=r_vals-1;
psi0=.5*(r_vals-r_0^2./r_vals)+sqrt(4*nu*t)*zeta_vals(:,1)-y_vals;
psi1=sqrt(4*nu*t)*zeta_vals(:,4);

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

% compute the streamfunction (note that this are column-row vector
% multiplications producing two-dimensional grids of psi-values.)
psi=psi0*u_e+psi1*(u_e.*u_ep)*(t/r_0);

% 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)))])
