%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Plot the two-term Blasius expansion for the vorticity field %%%%%
%%%%%%%%%% 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;

% 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 vorticity profiles
ome0=zeta_vals(:,3)/sqrt(4*nu*t);
ome1=zeta_vals(:,6)/sqrt(4*nu*t);

% compute the vorticity
ome=ome0*u_e+ome1*(u_e.*u_ep)*(t/r_0);

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

% plot the vorticity
contour(x,y,ome,-110:20:110)
