%	First approximation to flow past flat plate; error function solution
%	Compare with cubic solution in BSL and quadratic in Churchill
%	All solutions are self-similar
%	Dimensionless distance along plate is normalized with respect to
%	free stream velocity divided by kinematic viscosity to illustrate the
%  time or distance evolution of the boundary layer.
%	eta is equal to y/sqrt(nu*x/U)
clear all
clf reset 
eta=linspace(eps,5);
v1st=erf(eta./sqrt(8));	% 1st approximation;  error function, t=2x/U
eta8=eta./sqrt(8);
v2nd=erf(eta8-(1/sqrt(8)).*(1-0.5.*sqrt(pi).*erf(eta8)./eta8));	% 2nd approx. with convection
vbsl=1.5.*(eta./4.64)-0.5.*(eta./4.64).^3;	%	BSL cubic solution
for i=1:100,
   if eta(i)>4.64 
      vbsl(i)=1.0;
   end
end
vquad=2.*(eta./5.48)-(eta./5.48).^2;			%	quadratic solution
load blasiusd.dat;
plot(eta,v1st,'g-',eta,v2nd,'k-',eta,vbsl,'r:',eta,vquad,'b--'...
   ,blasiusd(:,1),blasiusd(:,2),'k*','LineWidth',3.0)
ylabel('Velocity')
xlabel('eta=y/sqrt(nu x/U)')
title('Boundary Layer Velocity Profiles')
text(1.5,0.30,'green (solid)= error function, w/o trans. convect.')
text(1.5,0.25,'black (solid)= error function, w/ trans. convect.')
text(1.5,0.20,'blue (dash) = quadratic')
text(1.5,0.15,'red (dot) = cubic')
text(1.5,0.10,'* = Blasius exact solution')
input('Press Enter or Return key to continue');
%	Make movie of boundary layer velocity profile evolution
%	t=2x/U
clear all
clf reset
nx=30;
x=linspace(0.001,1,nx);
y=linspace(0,5);
M=moviein(nx);
for i=1:nx,
   vx=erf(y./sqrt(8*x(i)));
   plot(y,vx,'k-')
   hold on
   ylabel('Velocity')
	xlabel('Distance Normal to Plate')
   title('Boundary Layer Velocity Profile Evolution')
   M(:,i)=getframe;
end
movie(M,3)
