function update_velocity(obj, t0)
	% update velocity and predicted position of active tracks at time index t0

	n_active		= obj.count();

	%% get available position history over last 3 timesteps
	% if not available, encode these as NaN
	re_pos_history	= nan(3, 3, n_active);
	trel			= -2 : 0;
	for it = 1 : 3
		if t0 + trel(it) >= 1
			re_pos_history(:,it,:)	= obj.get_active_property(1:n_active, t0+trel(it), 're_pos', true);
		end
	end

	%% get velocity based on available position history
	n_history		= sum(isfinite(re_pos_history(1,:,:)), 2);
	re_pos			= squeeze(re_pos_history(:,3,:));
	re_vel			= nan(3, n_active);
	re_pred			= nan(3, n_active);
	
	is_ok			= isfinite(re_pos_history(1,:,:));
	len3			= all(is_ok(1,1:3,:), 2);
	len2			= all(is_ok(1,2:3,:), 2) & ~len3;
	len1			= ~len2 & ~len3;

	% linearly extrapolate based on last three frames where possible
	re_vel3			= 1/2*(re_pos_history(:,3,:) - re_pos_history(:,1,:));
	re_pred3		= -2/3*re_pos_history(:,1,:) + 1/3*re_pos_history(:,2,:) + 4/3*re_pos_history(:,3,:);
	re_vel(:,len3)  = re_vel3(:, len3);
	re_pred(:, len3)= re_pred3(:,len3);

	% linearly extrapolate based on last two frames where necessary
	re_vel2			= re_pos_history(:,3,:) - re_pos_history(:,2,:);
	re_pred2		= re_pos_history(:,3,:) + re_vel2;
	re_vel(:, len2) = re_vel2(:,len2);
	re_pred(:, len2)= re_pred2(:,len2);

	%% interpolate from predictor field where there are no data available
	ok				= len2 | len3;
	if nnz(ok) > 3
		% interpolate from neighbours
		% build scatteredInterpolant
		re_pos1			= re_pos(:,len1).';
		vi_pred			= scatteredInterpolant(re_pos(1,ok)', re_pos(2,ok)', re_pos(3,ok)', re_vel(1,ok)', 'linear', 'linear');		
		% evaluate for each velocity component
		re_vel(1,len1)	= vi_pred(re_pos1);
		vi_pred.Values	= re_vel(2,ok).';
		re_vel(2,len1)	= vi_pred(re_pos1);
		vi_pred.Values	= re_vel(3,ok).';
		re_vel(3,len1)	= vi_pred(re_pos1);
		% predict position
		re_pred(:,len1) = re_pos(:,len1) + re_vel(:,len1);
	else
		% no history is available, assign zero velocity
		re_vel(:,len1)  = 0;
		re_pred(:,len1) = re_pos(:,len1);
	end

 	%% assign to tracks
% 	active			= obj.tracks(subset);
% 	re_vel_ary		= cat(2, active.re_vel);
% 	re_pred_ary		= cat(2, active.re_pred);
% 	re_vel_ary(:,:,t0) = re_vel;
% 	re_pred_ary(:,:,t0)= re_pred;
% 	re_vel_ary		= num2cell(re_vel_ary , [1 3]);
% 	re_pred_ary		= num2cell(re_pred_ary, [1 3]);
% 	[active.re_pred]= re_pred_ary{:};
% 	[active.re_vel ]= re_vel_ary{:};
% 	obj.tracks(subset) = active;

	new_active		= cell(1, n_active);
	for n = 1 : n_active
		trk = obj.active(n);
		trk.re_vel(:,:,t0)   = re_vel(:, n);
		trk.re_pred(:,:,t0+1)= re_pred(:, n);
		new_active{n} = trk;
	end
	obj.active		= cell2mat(new_active);
end