%% input parameters
clear;
clc;

% location of data
dset_id=3;
switch dset_id
case 1
	res_file	= 'L:\array\tcf\statistics\30C_CD.mat';
	line_width  = 1;
	marker		= 'x-';
	marker_corr = 'o-';
	fig_offset	= 0;
	cmap_lines	= brewermap(5, 'Set1');
case 2
	res_file	= 'L:\array\tcf\statistics\30C_ref.mat';
	line_width  = 2;
	marker		= '-';
	marker_corr = '-';
	fig_offset	= 100;
	cmap_lines	= 0.7*brewermap(5, 'Set1');
case 3
	res_file	= 'L:\array\tcf\statistics\30C_EF.mat';
	line_width  = 1;
	marker		= 'x-';
	marker_corr = 'o-';
	fig_offset	= 200;
	cmap_lines	= brewermap(5, 'Set1');
end

% TCF flow parameters
h			= 150;			% channel half height
n_y			= 256;			% number of bins in y direction
u_tau		= 8.432100;		% reference friction velocity
U_bulk		= 168.739875;   % reference bulk velocity, mm/s
nu			= 1;			% viscosity, mm^2/s
delta_nu	= nu / u_tau;
kappa		= 0.41;			% von Karman constant
B			= 5.0;

% plotting parameters
font_size	= 24;



%% load result file
res			= load(res_file);

%% plot mean velocity profile
figure(fig_offset+1);
clf;
for ic = 1 : 3
	plot(res.u_mean(:,ic)/U_bulk, res.y/h, marker, ...
		'linewidth', line_width, 'color', cmap_lines(ic,:));
	hold on;
end
pbaspect([3/4 1 1]);
ylabel('$y/h$', 'interpreter', 'latex', 'fontsize', font_size);
xlabel('$U / U_b$', 'interpreter', 'latex', 'fontsize', font_size);
set(gca, 'fontsize', font_size, 'ticklabelinterpreter', 'latex');
xlim([-0.05 1.2]);
set(gca,'xtick',-0.0:0.2:1.2);

%% plot velocity variance
figure(fig_offset+2);
clf;

% normal stresses
layout	= tiledlayout(1, 2);
nexttile;
for i = 1 : 3
	legstr	= sprintf('$\\langle u_%d'' u_%d'' \\rangle$', i, i); 
	plot(res.uu_prime(:,i,i) / u_tau^2, res.y/h, marker, ...
		'linewidth', line_width, ...
		'color', cmap_lines(i,:), ...
		'DisplayName', legstr);
	hold on;
	legstr	= sprintf('$\\langle u_%d'' u_%d'' \\rangle$ corrected', i, i); 
	plot(res.uu_prime_corr(:,i,i) / u_tau^2, res.y/h, marker_corr, ...
		'linewidth', line_width, ...
		'color', cmap_lines(i,:), ...
		'DisplayName', legstr);
end
pbaspect([3/4 1 1]);
ylabel('$y/h$', 'interpreter', 'latex', 'fontsize', font_size);
xlabel('$\langle u_i u_j \rangle^+$', 'interpreter', 'latex', 'fontsize', font_size);
set(gca, 'fontsize', font_size, 'ticklabelinterpreter', 'latex');
set(legend, 'location', 'east', 'fontsize', font_size, 'interpreter', 'latex');
xlim([0 10]);
title('Reynolds Normal Stresses', 'interpreter', 'latex', 'fontsize', font_size);

% reynolds shear stresses
nexttile;
for i = 1 : 3
	j		= mod(i,3)+1;
	legstr	= sprintf('$\\langle u_%d'' u_%d'' \\rangle$', i, j); 
	plot(res.uu_prime(:,i,j) / u_tau^2, res.y/h, marker, ...
		'linewidth', line_width, ...
		'color', cmap_lines(i,:), ...
		'DisplayName', legstr);
	hold on;
	legstr	= sprintf('$\\langle u_%d'' u_%d'' \\rangle$ corrected', i, j); 
	plot(res.uu_prime_corr(:,i,j) / u_tau^2, res.y/h, marker_corr, ...
		'linewidth', line_width, ...
		'color', cmap_lines(i,:), ...
		'DisplayName', legstr);
end
plot([-1 1], [-1 1], 'k-', 'linewidth', 2, 'DisplayName', '$1-y/h$');
pbaspect([3/4 1 1]);
ylabel('$y/h$', 'interpreter', 'latex', 'fontsize', font_size);
xlabel('$\langle u_i u_j \rangle^+$', 'interpreter', 'latex', 'fontsize', font_size);
set(gca, 'fontsize', font_size, 'ticklabelinterpreter', 'latex');
set(legend, 'location', 'east', 'fontsize', font_size, 'interpreter', 'latex');
xlim([-1 1]);
title('Reynolds Shear Stresses', 'interpreter', 'latex', 'fontsize', font_size);

%% plot velocity profile on semilog scale
figure(fig_offset+10);
clf;
rng		= res.y <= 0;

for ic = 1 : 3
	legstr	= sprintf('$\\langle U_%d \\rangle$', ic);
	semilogx((h+res.y(rng))/delta_nu, res.u_mean(rng,ic)/u_tau, marker, ...
		'linewidth', line_width, ...
		'color', cmap_lines(ic,:), ...
		'DisplayName', legstr);
	hold on;
end
pbaspect([4/3 1 1]);

% viscous sublayer
yp_vsl	= linspace(0,5);
plot(yp_vsl, yp_vsl, 'r-', 'linewidth', 2, 'DisplayName', 'Viscous Sublayer');

% log law
yp_log	= logspace(1, 3);
up_log	= 1/kappa * log(yp_log) + B;
legstr	= sprintf('$\\kappa = %0.2f, B = %0.2f$', kappa, B);
plot(yp_log, up_log, 'k--', 'linewidth', 2, 'DisplayName', legstr);

set(legend, 'location', 'northwest', 'fontsize', font_size, 'interpreter', 'latex');
xlabel('$(y+h)^+$', 'interpreter', 'latex', 'fontsize', font_size);
ylabel('$U^+$', 'interpreter', 'latex', 'fontsize', font_size);
set(gca, 'fontsize', font_size, 'ticklabelinterpreter', 'latex');
ylim([-0.1 25]);
xlim([0.5 h/delta_nu]);

%% plot reynolds normal stress profile on semilog scale
figure(fig_offset+11);
clf;
rng		= res.y <= 0;

for i = 1 : 3
	j		= i;
	legstr	= sprintf('$\\langle u_%d u_%d \\rangle$', i, j);
	semilogx((h+res.y(rng))/delta_nu, res.uu_prime(rng,i,j)/u_tau^2, marker, ...
		'linewidth', line_width, ...
		'color', cmap_lines(i,:), ...
		'DisplayName', legstr);
	hold on;
	legstr	= sprintf('$\\langle u_%d u_%d \\rangle$ corrected', i, j);
	semilogx((h+res.y(rng))/delta_nu, res.uu_prime_corr(rng,i,j)/u_tau^2, marker_corr, ...
		'linewidth', line_width, ...
		'color', cmap_lines(i,:), ...
		'DisplayName', legstr);
end
pbaspect([4/3 1 1]);

set(legend, 'location', 'northwest', 'fontsize', font_size, 'interpreter', 'latex');
xlabel('$(y+h)^+$', 'interpreter', 'latex', 'fontsize', font_size);
ylabel('$\langle u_i u_j\rangle^+$', 'interpreter', 'latex', 'fontsize', font_size);
set(gca, 'fontsize', font_size, 'ticklabelinterpreter', 'latex');
ylim([0 10]);
xlim([.5 h/delta_nu]);
