%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% track (or resume tracking) tracers through JHU turbulence database
addpath([pwd '/../../../turbmat']);

%% input parameters
dset_dir	= 'L:/shadow/hit/';
init_pat	= '%s/ic.mat';		% pattern of initial conditions file
dset_pat	= '%s/%d.mat';		% pattern of dataset files

%% load initial condition
init_file	= sprintf(init_pat, dset_dir);
fprintf('loading initial condition %s\n', init_file);
ic			= load(init_file);
auth_token  = ic.auth_token;
dset_token  = ic.dset_token;
s_interp	= ic.s_interp;
t_interp	= ic.t_interp;
g_interp	= ic.g_interp;
t			= ic.t;
x0			= ic.x0;
n_tracers	= size(x0,2);
n_t			= length(t);

fprintf('   tracking %d tracers over %d timesteps [%0.3f to %0.3f]\n', ...
		n_tracers, n_t, t([1 end]));

%% search for existing dataset
k			= 1;
while k < n_t
	dset_file	= sprintf(dset_pat, dset_dir, k);
	if ~exist(dset_file, 'file')
		break;
	end
	k			= k+1;
end
if k == 1
	error('no initial condition found');
end
if k == n_t
	error('tracking has completed');
end
k0			= k-1;
fprintf('starting tracking from step %d (t=%0.3f)\n', k0, t(k0+1));

%% initialise tracers
fprintf('load step %d\n', k0);
dset_file	= sprintf(dset_pat, dset_dir, k0);
ds			= load(dset_file);
x0			= ds.x;
assert(size(x0,2) == n_tracers);

%% track tracers over Lagrangian trajectory
for k = k0 : n_t-2
	%% track tracers to next timestep
	t0			= t(k+1);% k is a zero based index
	t1			= t(k+2);
	dt_track	= t1-t0;
	fprintf('[%4d/%4d] t=%6.3f\n', k, n_t-1, t0);
	fprintf('   get position...\n');
	x1			= getPosition(auth_token, dset_token, t0, t1, dt_track, s_interp, n_tracers, x0); 
	
	%% query data at next timestep
	fprintf('   get velocity...\n');
	u1			= getVelocity(auth_token, dset_token, t1, s_interp, t_interp, n_tracers, x1);
	fprintf('   get gradient...\n');
	A1			= getVelocityGradient(auth_token, dset_token, t1, g_interp, t_interp, n_tracers, x1);
	A1			= permute(reshape(A1, [3 3 n_tracers]), [2 1 3]);

	%% save to file
	fprintf('   save...\n');
	dset_file	= sprintf(dset_pat, dset_dir, k+1);
	ds			= struct('x', x1, 'u', u1, 'A', A1, 't', t1);
	save(dset_file, '-struct', 'ds');
	
	%% update
	x0			= x1;
	u0			= u1;
	A0			= A1;
	t0			= t1;
end