close all
clear all

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

jl0=29;
viscos=1/8000;

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

% read modelled turbulent kinetic energy (file 1)
load te1_zonal_pans.mat
te1=te1_zonal_pans;
te3d=reshape(te1,ni,nj,nk);

% read modelled turbulent dissipation (file 1)
load ed1_zonal_pans.mat
ed1=ed1_zonal_pans;
ed3d=reshape(ed1,ni,nj,nk);


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


% average in 1st (x) and 3rd (z) direction
umean=mean(mean(u3d,3),1);
temean=mean(mean(te3d,3),1);
edmean=mean(mean(ed3d,3),1);


% face coordinates
load yc_zonal_pans.dat
yc=yc_zonal_pans;

y_zonal_pans(1)=0; % node coordinates
for j=2:nj
% dy = cell width
   dy(j)=yc(j)-yc(j-1);
   y(j)=0.5*(yc(j)+yc(j-1));
end
dy(nj)=0;
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_zonal_pans.ps  -deps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(2)
% plot te
plot(yplus,temean)
plot(temean,yplus,'linewidth',2)
hold
axis([0 5 0 8000])
handle=gca;
set(handle,'fontsize',[20])
ylabel('y^+','fontsize',[20])
xlabel('k^+','fontsize',[20])

% zoom. 'Position',[left bottom width height])
% create smaller axes in top right, and plot on it
%axes('Position',[.1 .4 0.4 0.4])
axes('OuterPosition',[.49 .51 0.43 0.43])
plot(yplus,temean)
plot(temean,yplus,'linewidth',2)
hold
axis([0 5 0 40])
plot(temean,yplus,'o','linewidth',2)
handle=gca;
set(handle,'fontsize',[20])


print k_rans_zonal_pans.ps  -deps
