%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% track sinking spheres through JHU turbulence database
% using Runge-Kutta exponential integration scheme
% 
% following
% 
% Peter J. Ireland, T. Vaithianathan, Parvez S. Sukheswalla, Baidurja Ray, Lance R. Collins,
% Highly parallel particle-laden flow solver for turbulence research,
% Computers & Fluids, Volume 76, 2013
%
% see also
% http://turbulence.pha.jhu.edu/docs/Database-functions.pdf
function k1 = jhu_track_sinking(dset_dir)
addpath([pwd '/turbmat']);

%% input parameters
%dset_dir	= 'M:/jhu/inertial/w1.00e-01/';
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;
x1			= ic.x0;
n_tracers	= size(x1,2);
n_t			= length(t);
lambda		= ic.lambda;
w_min		= ic.w_min;		% min slip velocity due to gravitational settling
w_max		= ic.w_max;		% max slip 
h			= mean(diff(t)); % step size

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
dset_file	= sprintf(dset_pat, dset_dir, k0);
ds_prev		= load(dset_file);
k1			= k0;

%% integrate with runge-kutta scheme
while k1 < n_t
	%% query data at current timestep
	k2			= k1+1;
	fprintf('[%4d/%4d] t=%6.3f - %6.3f\n', k1, n_t, t(k1), t(k2));
	fprintf('   get velocity...\n');
	% initial velocity
	x1			= ds_prev.x;
	%v1			= ds_prev.v;
	u1			= ds_prev.u;
	v1			= u1		+ w_min(:);		% intial velocity estimate
	% intermediate predictor step
	h			= t(k2) - t(k1);
	x2_int		= x1 + v1 * h;	
	u2_int		= getVelocity(auth_token, dset_token, t(k2), s_interp, t_interp, n_tracers, x2_int);
	% timestepping of velocity with corrector step
	v2_int		= (u1 + u2_int)/2 + w_min(:);	% intermediate velocity
	x2			= x1 + v2_int*h/2;				% corrected particle position at time t2
	u2			= getVelocity(auth_token, dset_token, t(k2), s_interp, t_interp, n_tracers, x2);
	v2			= u2 + w_min(:);				% new particle velocity at time t2
	% velocity gradient
	fprintf('   get gradient...\n');
	A2			= getVelocityGradient(auth_token, dset_token, t(k2), g_interp, t_interp, n_tracers, x2);
	A2			= permute(reshape(A2, [3 3 n_tracers]), [2 1 3]);
	
	%% save to file
	fprintf('   save...\n');
	dset_file	= sprintf(dset_pat, dset_dir, k2);
	ds			= struct('x', x2, 'v', v2, 'u', u2, 'A', A2, 't', t(k2));
	save(dset_file, '-struct', 'ds');
	
	%% update
	ds_prev		= ds;
	k1			= k2;
	continue;
	%% debug
	figure(1);
	clf;
	plot3(x2(1,:), x2(2,:), x2(3,:), '.');
	axis equal tight xy;
	drawnow;
end
