function sigma=middle_state(vl,vr,tol);
%
% For Riemann initial data, vl, vr in Riemann coordinates.
% Gives the middle state for stiff 2x2 systems of conservation laws.
% 'tol': accuracy.
% The routine uses Newton iteration to find the middle state.
%
if nargin<3,
	tol=1e-6;
end;
err=tol+1;
sigma0=[vr(1)-vl(1);vr(2)-vl(2)];
n=0;
while err>tol&n<1000,
	A=dJack(sigma0);
	b=Jack(sigma0)+vl-vr;
	e=A\b;
	sigma1=sigma0-e;
	err=sum(abs(e));
	sigma0=sigma1;
	n=n+1;
end;
if n==1000,
	'Not converged middle_state after 1000 iterations'
%else
%	[b' n]
end;
sigma=sigma1;
