close all
clear all

ni=66;
nj=49;
nk=66;

% read u1 & transform u1 to a 3D array
load u1_hybrid.mat
u1=u1_hybrid;
u3d=reshape(u1,ni,nj,nk);

% x coordinate direction = index i, first index
% y coordinate direction = index j, second index
% z coordinate direction = index k, third index

% compute average u by averaging in x (i) and z (k) directions
umean(1:nj)=0;
for k=1:nk
for j=1:nj
for i=1:ni
   umean(j)=umean(j)+u3d(i,j,k)/ni/nk;
end
end
end


% face coordinates
load yc_hybrid.dat
yc=yc_hybrid;

y_hybrid(1)=0; % node coordinates
for j=2:nj-1
% dy = cell width
   dy(j)=yc(j)-yc(j-1);
   y_hybrid(j)=0.5*(yc_hybrid(j)+yc_hybrid(j-1));
end
dy(nj)=0;
y_hybrid(nj)=yc_hybrid(nj-1);
y=y_hybrid;
yplus=y*8000;

figure(1)
% plot u log
plot(yplus,umean)
hold
B=5.2;
kappa=0.4;
% find y+=30
for j=2:nj
   ulog(j)=log(yplus(j))/kappa+B;
   if yplus(j) < 30
      if yplus(j+1) > 30
         jj=j;
      end
   end
end
ulog(nj)=0;
plot(yplus(jj:nj-1),ulog(jj:nj-1),'+')
axis([1 8000 0 31])
handle=gca;
set(handle,'fontsize',[20])
set(handle,'xscale','log')
xlabel('y^+','fontsize',[20])
ylabel('U^+','fontsize',[20])
print u_log_hybrid.ps  -deps
