close all
clear all


dx=0.1;
dz=0.05;

ni=34;
nj=49;
nk=34;

viscos=1/5200;

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

load v1_IDD_PANS.mat
v1=v1_IDD_PANS;
v3d1=reshape(v1,ni,nj,nk); % v_2 velocity

load w1_IDD_PANS.mat
w1=w1_IDD_PANS;
w3d1=reshape(w1,ni,nj,nk); % v_3 velocity

load te1_IDD_PANS.mat
te1=te1_IDD_PANS;
te3d1=reshape(te1,ni,nj,nk); % modeled turbulent kinetic energy

load eps1_IDD_PANS.mat
eps1=eps1_IDD_PANS;
eps3d1=reshape(eps1,ni,nj,nk); % modeled turbulent dissipation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  read v_2 & transform v_2 to a 3D array (file 2)
load u2_IDD_PANS.mat
u2=u2_IDD_PANS;
u3d2=reshape(u2,ni,nj,nk);

load v2_IDD_PANS.mat
v2=v2_IDD_PANS;
v3d2=reshape(v2,ni,nj,nk);

load w2_IDD_PANS.mat
w2=w2_IDD_PANS;
w3d2=reshape(w2,ni,nj,nk);

load te2_IDD_PANS.mat
te2=te2_IDD_PANS;
te3d2=reshape(te2,ni,nj,nk);

load eps2_IDD_PANS.mat
eps2=eps2_IDD_PANS;
eps3d2=reshape(eps2,ni,nj,nk);



% merge files
u3d = vertcat(u3d1,u3d2);
v3d = vertcat(v3d1,v3d2);
w3d = vertcat(w3d1,w3d2);
te3d = vertcat(te3d1,te3d2);
eps3d = vertcat(eps3d1,eps3d2);


%8 files (34x49x34)  of instantaneous data have been merged;  the new size of the arrays is 8*34x49x34
ni=size(u3d,1);

% 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);
vmean=mean(mean(v3d,3),1);
wmean=mean(mean(w3d,3),1);
temean=mean(mean(te3d,3),1);
epsmean=mean(mean(eps3d,3),1);


% face coordinates
load yc.dat

y_zonal_sdes(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;

tauw=viscos*umean(2)/y(2);
ustar=tauw^0.5;
yplus=y*ustar/viscos;;

load LM_Channel_5200_mean_prof.dat
y_DNS=LM_Channel_5200_mean_prof(:,1);
yplus_DNS=LM_Channel_5200_mean_prof(:,2);
u_DNS=LM_Channel_5200_mean_prof(:,3);

load LM_Channel_5200_vel_fluc_prof.dat
u2_DNS=LM_Channel_5200_vel_fluc_prof(:,3);
v2_DNS=LM_Channel_5200_vel_fluc_prof(:,4);
w2_DNS=LM_Channel_5200_vel_fluc_prof(:,5);
uv_DNS=LM_Channel_5200_vel_fluc_prof(:,6);

k_DNS=0.5*(u2_DNS+v2_DNS+w2_DNS);

% find equi.distant DNS cells in log-scale
xx=0;
for i=1:40
   ii = find(yplus_DNS < 10^xx);
   i1=ii(end);
   jDNS(i)=i1;
   xx=xx+0.2;
end


% T.2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
figure(1)
% plot u log
plot(yplus,umean/ustar,'linewidth',2)
hold
h1=plot(yplus_DNS(jDNS),u_DNS(jDNS),'bo','linewidth',2)
set(h1,'MarkerFaceColor','white')
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.ps  -depsc

