%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute mean velocity profile statistics for TCF data
%

%% input parameters
%src_dir		= 'L:\array\tcf\30C\tracks\';
src_dir		= 'L:\array\tcf\rescaled\';
src_pat		= '%s/z%03d_t%03d.mat';
src_rng_z	= 1 : 1;
src_rng_t	= 1 : 1;
src_subset  = 1 : 100000;

h			= 150;			% channel half height
n_y			= 64;			% 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;

res_pat		= '%s/test.mat';
font_size	= 24;
fig_offset  = 100;

%% create result structure
[iz,it]			= ndgrid(src_rng_z, src_rng_t);

res				= struct();
res.n_y			= n_y;
res.phi			= linspace(pi, 0, n_y+1);
res.y_edges		= h*cos(res.phi);
res.y			= 1/2 * (res.y_edges(1:end-1) + res.y_edges(2:end));
res.src_dir		= src_dir;
res.src_pat		= src_pat;
res.src_rng_z	= src_rng_z;
res.src_rng_t	= src_rng_t;
res.iz			= iz;
res.it			= it;

%% evaluate velocity profile
res.u_sum		= zeros(res.n_y, 3);
res.uu_sum		= zeros(res.n_y, 3, 3);
res.n_sum		= zeros(res.n_y, 1);

n_samp			= numel(iz);
for is = 1 : n_samp
	%% load file
	src_file= sprintf(src_pat, src_dir, iz(is), it(is));
	fprintf('[z%03d t%03d] load %s\n', iz(is), it(is), src_file);
	if exist(src_file, 'file') == 0
		warning('%s missing, skipping\n', src_file);
		continue;
	end
	ds		= load(src_file);
	
	%% is experimental?
	if isfield(ds, 'tracks')
		ds.x	= [tracks.x];
		ds.u	= [tracks.u];
	else
		% process subset of refrence data
		ds.x	= ds.x(:,src_subset);
		ds.u	= ds.u(:,src_subset);
		ds.n_part=size(ds.x,2);
	end
	
	%% quantize y position into bins
	[n_local,i_bin]	= histc(ds.x(2,:), res.y_edges);
	
	%% accumulate in array
	res.n_sum	= res.n_sum + n_local(1, 1:end-1).';
	ok		= i_bin > 0 & i_bin <= res.n_y;
	for i = 1 : 3 
		res.u_sum(:,i)			= res.u_sum(:,i)	 + accumarray(i_bin(1,ok).', ds.u(i,ok)', [res.n_y 1]);
		for j = 1 : 3
			uiuj				= ds.u(i,:) .* ds.u(j,:);
			res.uu_sum(:,i,j)	= res.uu_sum(:,i,j) + accumarray(i_bin(1,ok).', uiuj(1,ok).', [res.n_y 1]);
		end
	end
end
res.u_mean		= res.u_sum ./ res.n_sum;
res.uu_mean		= res.uu_sum ./ res.n_sum;
res.u_var		= res.uu_mean - res.u_mean .* reshape(res.u_mean, [res.n_y 1 3]);

%% plot mean velocity profile
figure(fig_offset+1);
clf;
for ic = 1 : 3
	plot(res.u_mean(:,ic)/U_bulk, res.y/h, '.-');
	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.u_var(:,i,i) / u_tau^2, res.y/h, '.-' ,'DisplayName', legstr);
	hold on;
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.u_var(:,i,j) / u_tau^2, res.y/h, '.-' ,'DisplayName', legstr);
	hold on;
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,  '.-', '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([.1 h/delta_nu]);

%% plot velocity 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.u_var(rng,i,j)/u_tau^2,  '.-', 'DisplayName', legstr);
	hold on;
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([.1 h/delta_nu]);
