%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Display the Blasius flat plate streamlines %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% flow velocity
U=1;

% length of the plate
l=1;

% Renolds number based on the length; leave this small enough that the
% streamline displacement effects can still be seen
Re_l=1000;

% viscosity
nu=U*l/Re_l;

% maximum eta value in the Blasius solution
eta_max=.8/sqrt(nu*l/U);

% scaled wall shear
tau_wall=0.33206;

% integrate the Blasius ODE to find f (in column 1)
[eta_vals, f_vals] = ode45(@blaseq, [0 eta_max], [0; 0; tau_wall]);

% To get more accurate potential flow streamlines, this program writes
% the Blasius solution in so-called optimal coordinates (x',y') that
% are related to the physical coordinates (x,y) by the relation:
%                x' + i y' = sqrt(x + i y)

% x'-values
xp_vals=sqrt(l*[0:.01:1]);

% y'-values
yp_vals=.5*sqrt(nu/U)*eta_vals;

% Blasius streamfunction in terms of x' and y'
psi=sqrt(nu*U)*f_vals(:,1)*xp_vals;

% compute the physical locations
[xp,yp]=meshgrid(xp_vals,yp_vals);
z=(xp+i*yp).^2;
x=real(z);
y=imag(z);

% plot the streamlines, with increased density near the wall 
contour(x,y,psi,[0 .005 .01 .02 .04:.02:ceil(max(max(psi)))])

% evaluate the displacement thickness
disp=1.72*sqrt(nu/U*x(1,:));

% evaluate a rough boundary layer thickness
del=4*sqrt(nu/U*x(1,:));

% add them to the plot in green and red
hold on
plot(x(1,:),del,'r',x(1,:),disp,'b',...
     x(1,:),disp+.025,'g',x(1,:),disp+.05,'g',x(1,:),disp+.1,'g')
hold off
