%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute mean velocity profile statistics for TCF data
%
% apply compensation of variances for simultaneous independent errors
clear;
clc;
addpath([pwd '/../../misc']);

%% input parameters
src_dir		= 'L:\array\tcf\';
src_pat_A	= '%s/30C_E/tracks/z%03d_t%03d.mat';
src_pat_B	= '%s/30C_F/tracks/z%03d_t%03d.mat';
% src_pat_A	= '%s/rescaled/z%03d_t%03d.mat';
% src_pat_B	= '%s/rescaled/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			= 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;

n_bins_du	= 100;
re_tol		= 0.5;			% tolerance for finding corresponding particle, mm
du_std_est	= [2; 3; 30];	% estimate of standard deviation of errors, mm/s

res_pat		= '%s/statistics/30C_EF.mat';
res_file	= sprintf(res_pat, src_dir);

%% create result structure
[src_iz,src_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_A	= src_pat_A;
res.src_pat_B   = src_pat_B;
res.src_rng_z	= src_rng_z;
res.src_rng_t	= src_rng_t;
res.src_iz		= src_iz;
res.src_it		= src_it;

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

res.du_edges	= du_std_est * linspace(-3, 3, n_bins_du+1);
res.du_bins		= 1/2*(res.du_edges(:,2:end) + res.du_edges(:,1:end-1));
res.du_hist		= zeros(3, n_bins_du+1);
res.du_pdf		= zeros(3, n_bins_du);

n_samp			= numel(src_iz);
for is = 1 : n_samp
	%% load data
	src_file_A	= sprintf(src_pat_A, src_dir, src_iz(is), src_it(is));
	src_file_B	= sprintf(src_pat_B, src_dir, src_iz(is), src_it(is));
	fprintf('[z%03d t%03d]\n', src_iz(is), src_it(is));
	fprintf('A: %s\n', src_file_A);
	res_A		= load(src_file_A);
	fprintf('B: %s\n', src_file_B);
	res_B		= load(src_file_B);
	
	% get reconstruction limits
	if isfield(res_B, 'setup')
		re_lim	= res_B.setup.re_bb_lim;
	else
		re_lim	= [res_B.x_lim; res_B.y_lim; res_B.z_lim].';
	end

	%% represent tracks as arrays
	if isfield(res_A, 'tracks')
		x_A		= [res_A.tracks.x];
		u_A		= [res_A.tracks.u];
	else
		% handle reference data
		x_A		= res_A.x(:, src_subset);
		u_A		= res_A.u(:, src_subset);
	end
	if isfield(res_B, 'tracks')
		x_B		= [res_B.tracks.x];
		u_B		= [res_B.tracks.u];
	else
		x_B		= res_B.x(:, src_subset);
		u_B		= res_B.u(:, src_subset);
	end
	n_A			= size(x_A,2);
	n_B			= size(x_B,2);
	
	%% find nearest neighbours
	[iAnearB, dist]= knnsearch(x_A.', x_B.', 'K', 1');
	b_close		= dist < re_tol;
	iB			= find(b_close);
	matches		= [iAnearB(b_close,1), iB(:)];
	n_AB		= size(matches, 1);
	
	fprintf('A   has %d particles\n', n_A);
	fprintf('B   has %d particles\n', n_B);
	fprintf('A&B has %d particles in common\n', n_AB);

	%% compare 
	x_pair_A	= x_A(:,matches(:,1));
	u_pair_A	= u_A(:,matches(:,1));
	x_pair_B	= x_B(:,matches(:,2));
	u_pair_B	= u_B(:,matches(:,2));

	x_AB		= 1/2 * (x_pair_A + x_pair_B);
	u_AB		= 1/2 * (u_pair_A + u_pair_B);
	du_AB		= 1/2 * (u_pair_A - u_pair_B);
	
	%% quantize y position into bins
	[n_local,i_bin]	= histc(x_AB(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 
		%% mean and variance
		res.u_sum(:,i)			= res.u_sum(:,i)	 + accumarray(i_bin(1,ok).', u_AB(i,ok)', [res.n_y 1]);
		for j = 1 : 3
			uiuj				=  u_AB(i,ok) .*  u_AB(j,ok);
			duiduj				= du_AB(i,ok) .* du_AB(j,ok);
			res.uu_sum(:,i,j)	= res.uu_sum(:,i,j)   + accumarray(i_bin(1,ok).',   uiuj(1,:).', [res.n_y 1]);
			res.dudu_sum(:,i,j)	= res.dudu_sum(:,i,j) + accumarray(i_bin(1,ok).', duiduj(1,:).', [res.n_y 1]);
		end
		
		%% histograms of error
		res.du_hist(i,:)	= res.du_hist(i,:) + histc(du_AB(i,:), res.du_edges(i,:));
	end
end

%% finish averaging
res.u_mean		= res.u_sum  ./ res.n_sum;
res.du_sum		= res.du_sum ./ res.n_sum;
res.uu_mean		= res.uu_sum ./ res.n_sum;
res.dudu_mean	= res.dudu_sum ./ res.n_sum;

% correction for independent errors
res.uu_mean_corr= res.uu_mean - res.dudu_mean;
res.uu_prime	= res.uu_mean - res.u_mean .* reshape(res.u_mean, [res.n_y 1 3]);
res.uu_prime_corr=res.uu_mean_corr - res.u_mean .* reshape(res.u_mean, [res.n_y 1 3]);

% standard error of the mean
res.u_sem		= sqrt(cat(2, res.uu_prime(:,1,1), res.uu_prime(:,2,2), res.uu_prime(:,3,3)) ./ res.n_sum); 

% evaluate error PDF
for i = 1 : 3
	res.du_pdf(i,:)	= res.du_hist(i,1:end-1) ./ midpoint(res.du_edges(i,:), res.du_hist(i,1:end-1), 2);
end

%% save to file
fprintf('saving result to %s\n', res_file);
save(res_file, '-struct', 'res');

