%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% numerically integrate the rotation of ellipsoidal particles
% subject to a velocity gradient along a predefined trajectory
function jhu_track_ellipsoid(lambda, St_r)
%% input parameters
dset_dir	= 'L:\shadow\ref\';
init_pat	= '%s/ic.mat';					% pattern of initial conditions file
dset_pat	= '%s/%d.mat';					% pattern of dataset files
part_pat	= '%s/../AR%0.2f/%d.mat';	% pattern of integrated result
font_size	= 24;
%lambda		= 2;							% aspect ratio parameter
%St_r		= 0.1;							% rotational stokes number

% turbulence parameters
tau_eta		= 0.0424;						% kolmogorov timescale, DNS units
nu			= 0.000185;						% kinematic viscosity, dns units
rho_f		= 1;							% fluid density

%% define particle dimensions
a			= lambda;
c			= 1;
V_p			= 4/3*pi*a*c*c;						% particle volume
a_eq		= (a*c*c)^(1/3);					% equivalent sphere radius
tau_r		= St_r*tau_eta;						% rotational relaxation time
rho_p		= 15/2 * tau_r*nu / a_eq^2 * rho_f;	% particle density
m_p			= rho_p * V_p;						% particle mass
I_p			= ellipsoid_inertia(m_p, a, c, c);	% moment of inertia tensor

params		= struct('a', a, 'b', c, 'c', c', ...
					 'rho_p', rho_p, 'rho_f', rho_f, ...
					 'nu', nu, 'mu', nu*rho_f, ...
					 'St_r', St_r, 'lambda', lambda, ...
					 'volume', V_p, 'mass', m_p, 'inertia', I_p);

%% load simulation conditions
init_file	= sprintf(init_pat, dset_dir);
ic			= load(init_file);
x0			= ic.x0;
t			= ic.t;
n_particles	= size(x0,2);
n_t			= length(ic.t);

fprintf('tracking %d particles over %d timesteps\n', n_particles, n_t);

%% create directory
part_file	= sprintf(part_pat, dset_dir, lambda, 0);
part_dir	= fileparts(part_file);
if ~exist(part_dir, 'dir')
	fprintf('creating %s\n', part_dir);
	mkdir(part_dir);
end

%% create initial conditions for integration
ds_file			= sprintf(dset_pat, dset_dir, 0);
fprintf('[%4d/%4d] load %s\n', 0, n_t-1, ds_file);
ds0				= load(ds_file);
% create random initial coordinate system
rng(1); % not-so-random random orientation
R0				= randrotmat(n_particles);
% initialise rotation rate in equilibrium (torque free) state
A0_body			= mtimesx(R0, 't', mtimesx(ds0.A, R0));
[T0,J0]			= stokes_torque(A0_body, a, c);
Omega0_body		= J0 \ -T0*0;
Omega0_lab		= mtimesx(R0, reshape(Omega0_body, [3 1 n_particles]));

% encode state
state0			= stokes_body_encode_state(R0, Omega0_body);

%% save initial conditions
ds0.Omega		= Omega0_lab;
ds0.R			= R0;
ds0.params		= params;
save(part_file, '-struct', 'ds0');

%% integrate trajectories
for k1 = 1 : n_t-1
	%% load next time step
	k0				= k1-1;
	fprintf('[%4d/%4d] load %s\n', k1, n_t-1, ds_file);
	tic;
	ds_file			= sprintf(dset_pat, dset_dir, k1);
	ds1				= load(ds_file);
	
	%% numerically integrate
	% velocity gradient tensor and time derivative
	% we compute the terms of a taylor series expansion about 
	%	t0.5 = (t0+t1)/2
	fprintf('            integrating %0.3f to %0.3f\n', ds0.t, ds1.t);
	dt				=  ds1.t - ds0.t;
	dA				= (ds1.A - ds0.A) / dt;
	A05				= (ds0.A + ds1.A)/2;
	t05				= (ds1.t + ds0.t)/2;
	
	% numerically integrate
	odefun			= @(t,state) stokes_body_ode(t-t05, state, A05, dA, params.inertia, params.a, params.c, params.mu);
	opts			= odeset('RelTol', 1E-6, 'AbsTol', 1E-6);
	[t_int,state]	= ode113(odefun, [ds0.t t05 ds1.t], state0(:), opts);
	% decode state after numerical integration
	state1			= reshape(state(end,:), [12 n_particles]);
	[R1,~,Omega1]	= stokes_body_decode_state(state1);
	
	%% update state and save
	ds1.Omega		= Omega1;	% angular velocity in lab frame
	ds1.R			= R1;		% orientation of body in lab frame
	ds1.params		= params;
	ds0				= ds1;
	state0			= state1;
	
	part_file		= sprintf(part_pat, dset_dir, lambda, k1);
	fprintf('            save %s\n', part_file);
	save(part_file, '-struct', 'ds1');
	toc;
end