function [u,s]=states(vl,vr,sigma,delta);
%
%	Finds the step function approximation to the Riemann solution,
%	with a middle state given by 'sigma' and accuracy by 'delta'.
%   The speed of the wave is calculated by the functions 'sspeed' and .
%   'rspeed'
%
if nargin<4,
	delta=0.05;
end;
u=[]; s=[];
um=0.5*([vl(1)+sigma(1)+phi(sigma(1)),vl(2)+phi(sigma(1))]'+...
        [vr(1)-phi(sigma(2)),vr(2)-sigma(2)-phi(sigma(2))]');
if sigma(1)==0,
	;
elseif sigma(1)<0,
	if sigma(2)~=0,
		u=um;
	end;
	s=sspeed(1,um,sigma(1));
else
	N=max(ceil(sigma(1)/delta),2);
	sig=linspace(0,sigma(1),N);
	if sigma(2)~=0,
		u=[vl(1)+sig(2:N);vl(2)*ones(1,N-1)];
	else
		u=[vl(1)+sig(2:N-1);vl(2)*ones(1,N-2)];
	end;
	s=rspeed(1,vl,sig);
end;
if sigma(2)==0,
	;
elseif sigma(2)<0,
	s=[s sspeed(2,um,sigma(2))];
else
	N=max(ceil(sigma(2)/delta),2);
	sig=linspace(0,sigma(2),N);
	u=[u [um(1)*ones(1,N-2);um(2)+sig(2:N-1)]];
	s=[s rspeed(2,um,sig)];
end;

	
