function xb = vcmline(x1, x2, x, y, lm, M, t)

% xb = line(x1, x2, x(VCs,:), y(VCs), lm(VCs), M, t);
% called by vcmbound
%
% It calculate the boundary point crossing with a line

K = 4;

[y1, t] = vcm(x1, x, y, lm, M, t);
y2 = vcm(x2, x, y, lm, M, t);
while y1 * y2 >= 0
	x2 = (5*x2 + x1) / 6;
	if norm(x2 - x1) < 1e-5
		disp('line not cross');
		break;
	end
	y2 = vcm(x2, x, y, lm, M, t);
end;

while 1,
	d = x2 - x1;
	if norm(d) < 1e-5
		break;
	end;
	
	xt = repmat(x1, K-1, 1) + (1:K-1)'/K * d;
	yt = vcm(xt, x, y, lm, M, t);
	if y1 < 0, 
		idx = find(yt > 0);
	else
		idx = find(yt < 0);
	end
	if isempty(idx),
		x1 = xt(K-1, :);
	elseif idx(1) > 1,
		x1 = xt(idx(1)-1, :);
		x2 = xt(idx(1), :);
	else
		x2 = xt(1, :);
	end
end;
xb = (x1 + x2) / 2;