function matches = link_position_velocity(x1, u1, t1, x2, u2, t2, recon)
	% link track end-points in position-velocity space
	% 
	% Input:
	%	(x1,u1,t1)	3xN,3xN,1xN arrays of end points (i.e. data missing t>t1)
	%   (x2,u2,t2)  3xN,3xN,1xN arrays of start points (i.e. data missing t<t2)
	%	recon		track reconnection options structure, with fields
	%	T_search		search interval between ends of tracks
	% 	n_truncate		number of data points to truncate before reconnecting
	%	C_threshold		1 x n_tau vector of cost threshold values C(tau)
	%	tau				1 x n_tau time lag for each value in C_threshold
	%
	% Output:
	%	matches		n_match x 2 array indicating unique correspondences
	%		(x1,u1,t1) -> (x2,u2,t2)
	%	matches(:,1)	are indices in (x1,u1,t1)
	%	matches(:,2)	are corresponding points in (x2,u2,t2)
	% 
	% candidates are matched by minimising the generalised distance
	% 
	% D = |dx|^2 + |du|^2*|dt|^2/12
	% where
	%	dx = x2 - x1
	%	du = u2 - u1
	%	dt = t2 - t1
	
	n_x1		= size(x1,2);
	n_x2		= size(x2,2);
	assert(size(x1,1) == 3 && ismatrix(x1));
	assert(size(x2,1) == 3 && ismatrix(x2));
	assert(size(u1,1) == 3 && ismatrix(u1) && size(u1,2) == n_x1);
	assert(size(u2,1) == 3 && ismatrix(u2) && size(u2,2) == n_x2);
	assert(isvector(t1) && size(t1, 2) == n_x1);
	assert(isvector(t2) && size(t2, 2) == n_x2);
	
	%% find all possible connections
	% indices i0 
	ind2		= rangesearch(t2, t1, recon.T_search);
	n_candidates= cellfun(@(x) numel(x), ind2);
	ind2		= cell2mat(ind2.');
	
	ind1		= repelem(1 : n_x1, n_candidates);
	% indices satisfy | t2(ind1) - t1(ind2) | < T_search 
	% eliminate candidates where deltat < 0
	dt			= t2(1,ind2) - t1(1,ind1);
	ok			= dt > 0;
	candidates  = [ind1(1,ok); ind2(1,ok)].';

	%% evaluate cost function 
	% D = |dx|^2 + |du|^2*|dt|^2/12
	% where
	%	dx = x2 - x1
	%	du = u2 - u1
	%	dt = t2 - t1
	dx			= x2(:, candidates(:,2)) - x1(:, candidates(:,1));
	du			= u2(:, candidates(:,2)) - u1(:, candidates(:,1));
	dt			= t2(1, candidates(:,2)) - t1(1, candidates(:,1));
	cost		= sum(dx.^2 + (du.*dt).^2/12, 1);
	% evaluate threshold cost vector
	C_threshold	= interp1(recon.tau, recon.C_threshold, dt, 'linear', 0);
	accept		= cost < C_threshold;
	% create cost matrix
	cost_matrix = sparse(candidates(accept,1), candidates(accept, 2), 1./cost(1,accept), n_x1, n_x2);

	
	%% matches maximise 1./cost
	% the cost matrix defines a bipartite graph representing potential
	% matches between points X1=(x1,u1,t1) and X2=(x2,u2,t2)
	% our task is to find a matching of this graph which selects the best
	% possible candidates
	%
	% https://en.wikipedia.org/wiki/Maximum_cardinality_matching
	% https://en.wikipedia.org/wiki/Matching_(graph_theory)
	%
	% we solve this using a simple, greedy algorithm which discards ties 
	% handle ties to ensure there can only be a 1:1 correspondence of A->B
	[maxv,best]		= max(cost_matrix,[],1);
	iB				= find(maxv~=0);
	iA				= best(maxv~=0); 

	% assemble in array
	matches			= [iA(:), iB(:)];
end