%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare simultaneous independent measurements A and B
% made from different perspectives

%% input parameters
src_dir		= 'L:\array\tcf\';
src_pat_A	= '%s/6C/tracks/z%03d_t%03d.mat';
src_pat_B	= '%s/6Cm/tracks/z%03d_t%03d.mat';
src_iz		= 4;
src_it		= 1:1;
% tolerance for finding neighbours
re_tol		= 0.5;	% mm 

n_bins_du	= 100;

%% load measurements
src_file_A	= sprintf(src_pat_A, src_dir, src_iz, src_it);
src_file_B	= sprintf(src_pat_B, src_dir, src_iz, src_it);

fprintf('reconstruction A: %s\n', src_file_A);
res_A		= load(src_file_A);
fprintf('reconstruction B: %s\n', src_file_B);
res_B		= load(src_file_B);
re_lim		= res_B.setup.re_bb_lim;

% represent as arrays 
x_A			= [res_A.tracks.x];
u_A			= [res_A.tracks.u];
n_A			= size(x_A,2);

x_B			= [res_B.tracks.x];
u_B			= [res_B.tracks.u];
n_B			= size(x_B,2);

fprintf('A has %d particles\n', n_A);
fprintf('B has %d particles\n', n_B);

%% 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_matches	= size(matches, 1);

%% 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_mean		= 1/2 * (x_pair_A + x_pair_B);
u_sigma		= 1/2 * (u_pair_A + u_pair_B);
u_delta		= 1/2 * (u_pair_A - u_pair_B);

%% evaluate statistics of velocity difference
du_std		= std(u_delta,[],2);
du_mean		= mean(u_delta, 2);
du_rms		= rms(u_delta, 2);

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

for i = 1 : 3
	du_hist(i,:)	= histc(u_delta(i,:), du_edges(i,:));
	du_pdf(i,:)		= du_hist(i,1:end-1) ./ midpoint(du_edges(i,:), du_hist(i,1:end-1), 2);
end

du_mean_sem= du_std / sqrt(n_matches);

%% report
fprintf('found %d particles within %0.3f mm\n', n_matches, re_tol);
fprintf('%4.1f%% of A\n', n_matches / n_A * 100);
fprintf('%4.1f%% of B\n', n_matches / n_B * 100);

fprintf('RMS velocity error\n');
fprintf('u    %6.2f mm/s\n', du_rms(1));
fprintf('v    %6.2f mm/s\n', du_rms(2));
fprintf('w    %6.2f mm/s\n', du_rms(3));
fprintf('mean velocity error\n');
fprintf('u    %6.2f mm/s +- %4.2f\n', du_mean(1), 2*du_mean_sem(1));
fprintf('v    %6.2f mm/s +- %4.2f\n', du_mean(2), 2*du_mean_sem(2));
fprintf('w    %6.2f mm/s +- %4.2f\n', du_mean(3), 2*du_mean_sem(3));

%% plot particle distribution
figure(1);
clf;
plot3(x_A(1,:), x_A(2,:), x_A(3,:), 'r.');
hold on;
plot3(x_B(1,:), x_B(2,:), x_B(3,:), 'b.');
plot3(x_mean(1,:), x_mean(2,:), x_mean(3,:), 'ko');

axis equal tight xy;
axis(re_lim(:));

%% plot histogram of velocity error
figure(2);
clf;
cmap	= eye(3);
for i = 1 : 3
	semilogy(du_bins(i,:) / du_std(i), du_pdf(i,:), '.-', 'color', cmap(i,:));
	hold on;
	%plot(du_mean(i)*[1 1], [1e-4 1], '--', 'color', cmap(i,:));
end
pbaspect([4/3 1 1]);
xlabel('$\Delta u / \sigma$', 'interpreter', 'latex', 'fontsize', font_size);
ylabel('$P(\Delta u)$', 'interpreter', 'latex', 'fontsize', font_size);
set(gca, 'ticklabelinterpreter', 'latex', 'fontsize', font_size);
