close all
clear all


%************
viscos=1.52e-5;
i_lower_sym=20; % the flat plate starts at i=i_lower_sym+1
load xc.dat
load yc.dat
load boundary_layer_data.dat
file=boundary_layer_data;
n=length(file);

ni=252;
nj=200
n=length(file);

u=file(:,1);  %don't use this array
v=file(:,2);  %don't use this array
p=file(:,3);  %don't use this array


u2d=reshape(u,nj,ni)'; % this is v_1
v2d=reshape(v,nj,ni)'; % this is v_2
p2d=reshape(p,nj,ni)'; % this is p


% scale u and v with max u
for i=1:ni
   u2d(i,:)=u2d(i,:)/max(u2d(i,:));
   v2d(i,:)=v2d(i,:)/max(u2d(i,:));
end


load blasius.dat
xi_blas=blasius(:,1);
g_blas=blasius(:,2);
u_blas=blasius(:,3);


%   a control volume, CV. 
%
%  xp(i), yp(j) denote the center of the, CV. u, v and p are stored at (xp,yp)
%
%  xc(i) yc(j) denote the corner (on the high side) of the CV
%
%
%   x-------------------------x  xc(i), yc(j)
%   |                         |
%   |                         |
%   |                         |
%   |                         |
%   |          o xp(i), yp(j) |
%   |                         |
%   |                         |
%   |                         |
%   |                         |
%   x-------------------------x
%
% compute xp
xp(1)=xc(1);
for i=2:ni-1
   xp(i)=0.5*(xc(i)+xc(i-1));
end
xp(ni)=xc(ni-1);

% compute yp
yp(1)=yc(1);
for j=2:nj-1
   yp(j)=0.5*(yc(j)+yc(j-1));
end
yp(nj)=yc(nj-1);
%
% make xp and yp 2D arrays
for i=1:ni
   y2d(i,:)=yp;
end
for j=1:nj
   x2d(:,j)=xp;
end



%
% compute the gradient dudx, dudy
[dudy dudx]=gradient(u2d,yp,xp);


%************
figure(1)

%************
% velocity profile plot
figure(1)
i=170; % plot the velocity profile for i=170
plot(u2d(i,:),yp,'linewidth',2);
hold  % hold the plot 
i=5; % plot the velocity profile for i=5
plot(u2d(i,:),yp,'r--','linewidth',2);  %red dashed line with line width 2
title('Velocity profile','fontsize',[20])
axis([0 1.1 0 0.05]) % set x & y axis
xlabel('V_1','fontsize',[20]) %20-pt text
ylabel('x_2','fontsize',[20]) %20-pt text
text(0.04,0.01,'x_1=0.008 & 0.34','fontsize',[20]) % show this text at (0.04,0.004)
handle=gca;
set(handle,'fontsize',[20]); %20-pt text on the axis
print -deps velprof.ps

%************
% contour plot of v1
figure(2)
surf(x2d,y2d,u2d);
shading interp  % choose type of shading
view(0,90)      % look at the plot from the positive z-axis (2D)
axis('equal')   % make the x and y axis equal, i.e plot the channel in correct proportion
axis([0 0.1 0 0.011]) % zoom-in on the first 0.1m from the inlet
caxis([0.1 2.]) % set color axis [min,max]
handle1=colorbar('SouthOutside'); %colorbar at the bottom
set(handle1,'Position',[0.13 0.3 0.775 0.008]);  % colorbar(x,y,width,height)
handle=gca;
set(handle,'fontsize',[20]); %20-pt text on the colorbar
title('V1','fontsize',[20])
print -depsc v1_iso.ps


