function [u,s]=Riemannsol(vl,vr,delta);
%
%  Finds the Riemann solution for a 2x2 system (isentropic p-system, gamma=1)
%  (vl,vr) initial data in RIEMANN COORDINATES, delta: size of steps in approximation.
%
if max(abs(vl-vr))<delta*0.0001, % this problem is just too small...
	s=[]; u=vl;
else
	sigma=middle_state(vl,vr,delta/1000.0); % Finds the middle state 
	[u s]=states(vl,vr,sigma,delta); % Finds the states in the two waves
end;
