
t_end=2.8;
N=70+1;
J=1000;

h = 2.0 / (N-1);
x = [-1.0:h:1.0];

%%%
% different initial values
%%%

% two shocks moving to the left
u01 =-0.4*ones(size(x)) - 0.9.*( x > 0.5 ) + 0.2 *( x <-0.2 );

% one shock and one rarefaction moving to the right
u02 = 0.2*ones(size(x)) + 0.9.*( x <-0.5 ) + 0.2 *( x > 0.2 );

% shock moving to the right and transsonic rarefaction
u03 =-0.3*ones(size(x)) + 1.3.*( x <-0.5 ) + 0.4 *( x > 0.2 );

% transsonic rarefaction and shock moving to the right
u04 = 0.2*ones(size(x)) - 0.5.*( x <-0.3 ) - 0.2 *( x > 0.5 );

% single transsonic rarefaction
u05 = 0.3*ones(size(x)) - (x<0. ) - 0.5*(x==0. );

U = u04;


figure(1)
plot(x,U,'o-')
u_min = min(U)-.2;
u_max = max(U)+.2;
xlim([-1 1]); ylim([u_min u_max]);
xlabel('x'); ylabel('u');
pause

% u_t + a u_x = 0
% shock-only Riemann solver
% with  entropy fix for trans-sonic rarefaction
t = .0;
for ind = 1:J
	f = 0.5*U.^2;
	df = U;
	ur_ul = U(2:end) -U(1:end-1);
  	ur_ul(ur_ul==.0) = 1.0;
	s = (f(2:end) -f(1:end-1))./ur_ul  +df(2:end).*(ur_ul==.0);
	W = U(2:end) -U(1:end-1);
	tau = .75*h/max(abs(s));
	F_r = [min(s,.0).*W 0];
	F_l =-[0 max(s,.0).*W];

	% entropy fix
	idx = find( df(1:end-1)<.0 & df(2:end)>.0 );
	F_r(idx  ) =-f(idx);
	F_l(idx+1) =-f(idx+1);

	U = U -tau/h* ( F_r -F_l ) ;

	plot(x,U,'o-')
	xlim([-1 1]); ylim([u_min u_max]);
	xlabel('x'); ylabel('u');
	drawnow

	t = t + tau;
	if(t >= t_end) break; end;
end
