% Channel flow
close all
clear all
load channel_flow_data.dat

ni=199;  % number of grid nodes in x_1 direction
nj=28;  % number of grid nodes in x_2 direction
x1=channel_flow_data(:,1); %don't use this array
x2=channel_flow_data(:,2); %don't use this array
v1=channel_flow_data(:,3); %don't use this array
v2=channel_flow_data(:,4); %don't use this array
p=channel_flow_data(:,5);  %don't use this array

% transform the arrays from 1D fields x(n) to 2D fields x(i,j)
% the first index 'i', correponds to the x-direction
% the second index 'j', correponds to the y-direction

x1_2d=reshape(x1,ni,nj); %this is x_1 (streamwise coordinate)
x2_2d=reshape(x2,ni,nj); %this is x_2 (wall-normal coordinate)
v1_2d=reshape(v1,ni,nj); %this is v_1 (streamwise velocity component)
v2_2d=reshape(v2,ni,nj); %this is v_1 (wall-normal velocity component)
p_2d=reshape(p,ni,nj);   %this is p   (pressure)


%************
% velocity profile plot
figure(1)
i=170; % plot the velocity profile for i=170
plot(v1_2d(i,:),x2_2d(i,:),'linewidth',2);
hold  % hold the plot 
i=5; % plot the velocity profile for i=5
plot(v1_2d(i,:),x2_2d(i,:),'r--','linewidth',2);  %red dashed line with line width 2
title('Velocity profile','fontsize',[20])
axis([0 1.5 0 0.01]) % set x & y axis
xlabel('V_1','fontsize',[20]) %20-pt text
ylabel('x_2','fontsize',[20]) %20-pt text
text(0.04,0.004,'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(x1_2d,x2_2d,v1_2d);
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

%************
figure(3)
% compute the velocity gradient dv_1/dx_2
for i=1:ni
for j=2:nj-1
   dx2=x2_2d(i,j+1)-x2_2d(i,j-1);
   dv1_dx2(i,j)=(v1_2d(i,j+1)-v1_2d(i,j-1))/dx2;
end
end

% fix the derivative at the walls
for i=1:ni
% lower wall
   dx2=x2_2d(i,2)-x2_2d(i,1);
   dv1_dx2(i,1)=(v1_2d(i,2)-v1_2d(i,1))/dx2;
% upper wall
   dx2=x2_2d(i,nj)-x2_2d(i,nj-1);
   dv1_dx2(i,nj)=(v1_2d(i,nj)-v1_2d(i,nj-1))/dx2;
end

i=170; % plot the velocity gradient for i=170
plot(dv1_dx2(i,:),x2_2d(i,:),'linewidth',2);
hold
axis([-550 550 0 0.01]) % set x & y axis
title('Velocity gradient','fontsize',[20])
xlabel('dv_1/dx_2','fontsize',[20])
ylabel('x_2','fontsize',[20]) 
text(-380,0.004,'x_1=0.52','fontsize',[20])
handle=gca;
set(handle,'fontsize',[20]);
print -deps dv1_dx2.ps

