%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Plot the two-term Blasius expansion for the velocity profile %%%%
%%%%%%%%%% 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

% values of boundary layer coordinate y
r_max=0.15*r_0;
y=[0:.003*r_0:r_max];

% angular position
thet=30
% u_e
u_e=2*U*sin(thet*pi/180);
% its x-derivative
u_ep=2*U*cos(thet*pi/180)/r_0;
% the velocity values according to Blasius
u1=u_e*erf(y./sqrt(4*nu*t))+t*u_e*u_ep*zeta1p(y./sqrt(4*nu*t))...
   +.5*u_e*(1+r_0^2./(r_0+y).^2)-u_e;

% angular position
thet=90
% u_e
u_e=2*U*sin(thet*pi/180);
% its x-derivative
u_ep=2*U*cos(thet*pi/180)/r_0;
% the velocity values according to Blasius
u2=u_e*erf(y./sqrt(4*nu*t))+t*u_e*u_ep*zeta1p(y./sqrt(4*nu*t))...
   +.5*u_e*(1+r_0^2./(r_0+y).^2)-u_e;

% angular position
thet=150
% u_e
u_e=2*U*sin(thet*pi/180);
% its x-derivative
u_ep=2*U*cos(thet*pi/180)/r_0;
% the velocity values according to Blasius
u3=u_e*erf(y./sqrt(4*nu*t))+t*u_e*u_ep*zeta1p(y./sqrt(4*nu*t))...
   +.5*u_e*(1+r_0^2./(r_0+y).^2)-u_e;

% plot all three profiles in a single plot
u=[u1;u2;u3];
plot(u,y)
