%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialise tracers for tracking through JHU turbulence database
addpath([pwd '/turbmat']);

%% dimensional parameters of TCF dataset
delta_nu	= 1.0006E-3;		% plus unit
Ub			= 0.99994;			% bulk velocity
Uc			= 1.1312;			% centerline velocity
u_tau		= 4.9968E-2;		% friction velocity
nu			= 5E-5;				% kinematic viscosity
T_sim		= [0 25.9935];		% simulation duration
t_sim_step	= 0.0013;			% simulation timestep
t_save_step	= 0.0065;			% timestep for snapshots
h			= 1;				% half height
load('L:/array/tcf/grid.mat');

%% input parameters
dset_dir	= 'L:/array/tcf/';
dset_pat	= '%s/z%03d_t%03d.mat';	% pattern of initial conditions file

dt			= t_save_step*40;		% timestep for snapshots
t0			= 0;					% start time
t1			= 25.9935;				% end time

x_lim		= [0 8*pi];
y_lim		= [-1 1];
n_z			= 32;

% rough estimate for 0.005 part/px seeding concentration
FOV_h		= 2;
FOV_w		= 2 * 4/3;
FOV_d		= 50/300 * FOV_h;
rho_p		= 0.005 * 12e6 / (FOV_h * FOV_w * FOV_d);

% JHU access
auth_token	= 'uk.ac.cam.dm529-ca4ff15b';
dset_token  = 'channel';		% dataset
s_interp	= 'Lag4';			% 4th order Lagrangian interpolation in space
g_interp	= 'Fd4Lag4';		% interpolation for GetVelocityGradient
t_interp	= 'None';			% piecewise cubic interpolation in time

%% create time axis
t			= t0 : dt : t1;
n_t			= length(t);

% number of particles to query
FOV_x		= x_lim*[-1;1];
FOV_y		= y_lim*[-1;1];
FOV_z		= 3*pi/32;
n_part		= round(rho_p * FOV_x*FOV_y*FOV_z);
FOV			= [FOV_x FOV_y FOV_z]';

fprintf('Slab dimensions: %0.3f x %0.3f x %0.3f\n', FOV);
fprintf('rho_p:           %g\n', rho_p);
fprintf('# particles:     %d\n', n_part);

%% iterate over timesteps and slab positions
for iz = 1 : n_z
	z_offset = 3*pi / n_z * (iz-1);
	z_lim	= z_offset + [0 3*pi/n_z];
	
	for it = 1 : n_t
		%% create distribution of tracers
		x0			= diag(FOV) * rand(3,n_part) + [x_lim(1); y_lim(1); z_lim(1)];

		%% create initialisation files
		fname		= sprintf(dset_pat, dset_dir, iz, it);
		fs			= struct();
		fs.t		= t(it);
		fs.x		= x0;
		fs.x_lim	= x_lim;
		fs.y_lim	= y_lim;
		fs.z_lim	= z_lim;
		fs.n_part   = n_part;
		fs.s_interp	= s_interp;
		fs.g_interp	= g_interp;
		fs.t_interp = t_interp;
		fs.auth_token=auth_token;
		fs.dset_token=dset_token;
		fprintf('[z%03d, t%03d] z = %0.3f-%0.3f, t = %0.3f\n', iz, it, z_lim, t(it));
		fprintf('%s\n', fname);
		save(fname, '-struct', 'fs');
	end
end
