function [u,x]=genFrontTrack(u0,x0,T,delta,varargin);
%
%	Does front tracking solution of the initial data 'u0(1:n)', 
%   with discontinuities at the points 'x0(1:n-1)'up to time T, 
%   using 'delta' as a parameter in the Riemann solver.
%
%   The following parameters can be set: 'display':'lines'      : for drawing the fronts as lines
%										    	   'polygons'   : for drawing colored polygons
%                                                  'none'      : no graphics, default.
%
%                                       'save'    : creates output that can be saved, default not.
%
% 										'colors'  : a colormap to be used for colors, must be
%												    followed by the colormap.
%
%                                       'periodic': if periodic boundary, must be followed by
%                                                   the interval '[a,b]', where a<x0<b. Also
%                                                   u0(:,1)=u0(:,length(u0)) in this case.
%  
%	IMPORTANT:
%   To use this you must supply our own routine [u s]=Riemannsol(ul,ur,delta);


%% Preliminaries ....
l=length(varargin);
plain=0; lines=0; polygons=0; savegraphics=0; cmap_set=0; Periodic=0; COLL_TOL=0.0;
if l==0,
	plain=1;
else
	i=1;
	while i<=l
		if strcmp(varargin(i),'display')
			i=i+1;
			if strcmp(varargin(i),'polygons')
				polygons=1;
			elseif strcmp(varargin(i),'lines')
				lines=1;
			elseif strcmp(varargin(i),'none')
				plain=1;
			else
				error('                  No such option for display');
			end;
		elseif strcmp(varargin(i),'save')
			savegraphics=1;
		elseif strcmp(varargin(i),'colors')
			i=i+1;
			cmap=varargin(i);
			cmap=cmap{1,1};
			cmap_set=1;
		elseif strcmp(varargin(i),'periodic')
			Periodic=1;
			i=i+1;
			bndr=varargin(i);
			xleft=bndr{1}(1);
			xright=bndr{1}(2);
		else
			error('                    No such code')
		end;
		i=i+1;
	end;
end;
%%
% To avoid lots of 'if' tests I chose to create rather bulky code
if plain,     % READ THIS CASE IF CONCERNED WITH FRONT TRACKING, the other cases differ only in 
	% graphics option, such as plotting in various modes. 
	% Initializes the fronts and collision times
	if ~Periodic,
		[u s t_coll,x]=initialize(u0,x0,T,delta);
	else
		[u s t_coll,x]=initialize(u0,x0,T,delta,'periodic',xleft,xright);
	end;
	% 'u' are the states AFTER the initial Riemann problems have been solved
	% 's' are the speeds of the discontiuity between 'u(i-1)' and 'u(i)'
	% 't_coll' is the collision time of fronts at 'x(i-1)' and 'x(i)'
	t_start=zeros(size(x));
	% The starting times of the fronts
	if Periodic,
		xlower=xleft;
		xhigher=xright;
	else
		xlower=min([x(1)+s(1)*T,x(1)]);
		xhigher=max([x(length(x))+s(length(s))*T,x(length(x))]);
		if xlower==xhigher
			xlower=xhigher-0.5;
			xhigher=xlower+1;
		end;
	end;
	[t i]=min(t_coll);        
	% The first collision time is between fronts i-1 and i
	first=1;
	% Continues tracking while there are more collsions before T, and more
	% discontinuities
	ncoll=0;
	while t<T&~isempty(s), % Tracking 
		ncoll=ncoll+1;
		[ncr nu]=size(u);
		nx=length(x);
 		if nx~=nu-1,
			[nx nu];
			error('Wrong length of x vs u in FrontTrack');
		end;
		% Finds the collsion point
		xc=x(i)+s(i)*(t-t_start(i));
		if i<nu-2,
			xp=[x(i) x(i+1) xc];
			tp=[t_start(i) t_start(i+1) t];
		else
			xp=[x(i) xhigher xhigher xc];
			tp=[t_start(i) t_start(i) t t];
		end;
		% solving the Riemann problem, giving intermediate states 'ur' 
		% and speeds 'sr'
		if (Periodic&i==2),
			tcoll_left=collision_time([0.0 s(3)], [t t_start(3)], [xleft x(3)],T+1);
			tcoll_right=collision_time([s(nu-2) s(2)], [t_start(nu-2) t],...
				[x(nu-2) xright],T+1);
			x=[xleft x(3:nu-2) xright xright];
			s=[0 s(3:nu-2) s(2) 0];
			t_start=[t t_start(3:nu-2) t t];
			u=[u(:,3) u(:,3:nu-1) u(:,3) u(:,3)];
			t_coll=[T tcoll_left t_coll(4:nu-2) tcoll_right T T];
		elseif (Periodic&i==nu-1),
			tcoll_right=collision_time([s(nx-2) 0.0],[t_start(nx-2) t],...
				[x(nx-2) xright],T+1);
			tcoll_left=collision_time([s(nx-1),s(2)],[t t_start(2)],...
				[xleft, x(2)],T+1);
			x=[xleft xleft x(2:nx-2) xright];
			t_start=[t t t_start(2:nx-2) t];
			s=[0 s(nx-1) s(2:nx-2) 0];
			u=[u(:,nu-2) u(:,nu-2) u(:,2:nu-2) u(:,nu-2)];
			t_coll=[T,T, tcoll_left,t_coll(3:nu-3),tcoll_right,T];
		else
			[ur sr]=Riemannsol(u(:,i-1),u(:,i+1),delta);  
			[nc nr]=size(ur); 
			ns=length(sr);
			% updates the collision vector
			if ns>0, 
				if i>2,
					t_coll(i-1)=collision_time([s(i-2) sr(1)],[t_start(i-2),t],...
					[x(i-2),xc],T+1);
					if t_coll(i-1)<t,
						[xc t_coll(i-1) t]
						[s(i-2) sr(1) t_start(i-2) t]
						error('not increasing t_coll');
					end;
				end;
				if i<nu-1,
					t_coll(i+1)=collision_time([sr(ns),s(i+1)],[t,t_start(i+1)],...
					[xc,x(i+1)],T+1);
					if t_coll(i+1)<t,
						[xc t_coll(i+1) t]
						[sr(ns) s(i+1)]
						error('not increasing t_coll');
					end;
				end;
			else
				ur=[]; sr=[]; nr=0;
				if i>2&i<nu-1,
					t_coll(i-1)=collision_time([s(i-2) s(i+1)],...
					[t_start(i-2),t_start(i+1)],...
					[x(i-2),x(i+1)],T+1);
				end;					  
			end;
			% Updates the speeds
			s=[s(1:i-2) sr s(i+1:nu-1)];   
			hone=ones([1,ns]);
			% Updates the locations
			x=[x(1:i-2) xc*hone x(i+1:nu-1)]; 
			% Updates the starting points
			t_start=[t_start(1:i-2) t*hone t_start(i+1:nu-1)]; 
			% Updates the values and the collsion times
			if ns,
				u=[u(:,1:i-1) ur u(:,i+1:nu)];
				t_coll=[t_coll(1:i-1) (T+1)*ones([1,nr]) t_coll(i+1:nu)];
			else
				u=[u(:,1:i-1) u(:,i+2:nu)];
				t_coll=[t_coll(1:i-1) t_coll(i+2:nu)];
			end;
		end;
		% Finds the next collision
		[t i]=min(t_coll);   
	end;
	% Moves the discontinuities up to T
	x=x+s.*(T-t_start); 	
else
	%% Drawing with simple lines
	if lines,
		if savegraphics, 
			%% Saving the graphics, takes some time and memory...
			scn=get(0,'ScreenSize');
			ext=[0.2 0.7];
			pos=[scn(1:2)+ext(1)*(scn(3:4)-scn(1:2)) scn(1:2)+ext(2)*(scn(3:4)-scn(1:2))];
			figure('Position',pos);
			if ~Periodic,
				[u s t_coll,x]=initialize(u0,x0,T,delta);
			else
				[u s t_coll,x]=initialize(u0,x0,T,delta,'periodic',xleft,xright);
			end;		
			t_start=zeros(size(x));
			if Periodic,
				xlower=xleft;
				xhigher=xright;
			else
				xlower=min([x(1)+s(1)*T,x(1)]);
				xhigher=max([x(length(x))+s(length(s))*T,x(length(x))]);
				if xlower==xhigher
					xlower=xhigher-0.5;
					xhigher=xlower+1;
				end;
			end;		
			axes('XLim',[xlower xhigher],'YLim',[0 T],'Box','on',...
				'Drawmode','fast');
			hold on;
			[t i]=min(t_coll);
			while t<T&~isempty(s),
				[ncr nu]=size(u);
				nx=length(x);
				if nx~=nu-1,
					[nx nu];
					error('Wrong length of x vs u in FrontTrack');
				end;
				xc=x(i)+s(i)*(t-t_start(i));
				line([x(i) xc],[t_start(i) t]);
				line([x(i-1) xc],[t_start(i-1) t]); 
				drawnow;
				if (Periodic&i==2),
					tcoll_left=collision_time([0.0 s(3)], [t t_start(3)], [xleft x(3)],T+1);
					tcoll_right=collision_time([s(nu-2) s(2)], [t_start(nu-2) t],...
						[x(nu-2) xright],T+1);
					x=[xleft x(3:nu-2) xright xright];
					s=[0 s(3:nu-2) s(2) 0];
					t_start=[t t_start(3:nu-2) t t];
					u=[u(:,3) u(:,3:nu-1) u(:,3) u(:,3)];
					t_coll=[T tcoll_left t_coll(4:nu-2) tcoll_right T T];
				elseif (Periodic&i==nu-1),
					tcoll_right=collision_time([s(nx-2) 0.0],[t_start(nx-2) t],...
						[x(nx-2) xright],T+1);
					tcoll_left=collision_time([s(nx-1),s(2)],[t t_start(2)],...
						[xleft, x(2)],T+1);
					x=[xleft xleft x(2:nx-2) xright];
					t_start=[t t t_start(2:nx-2) t];
					s=[0 s(nx-1) s(2:nx-2) 0];
					u=[u(:,nu-1) u(:,nu-1) u(:,2:nu-2) u(:,nu-2)];
					t_coll=[T,T, tcoll_left,t_coll(3:nu-3),tcoll_right,T];
				else
					[ur sr]=Riemannsol(u(:,i-1),u(:,i+1),delta);
					[nc nr]=size(ur); 
					ns=length(sr);
					if ns>0,
						if i>2,
							t_coll(i-1)=collision_time([s(i-2) sr(1)],[t_start(i-2),t],...
							[x(i-2),xc],T+1);
							if t_coll(i-1)<=t,
								[-1 t_coll(i-1) t]
								[s(i-2) sr(1) t_start(i-2) t]
								error('not increasing t_coll');
							end;
						end;
						if i<nu-1,
							t_coll(i+1)=collision_time([sr(ns),s(i+1)],[t,t_start(i+1)],...
							[xc,x(i+1)],T+1);
							if t_coll(i+1)<=t,
								[1 t_coll(i+1) t]
								[sr(ns) s(i+1)]
								line([xc,xc+sr(ns)*(T-t)],[t,T]);
								drawnow;
							end;
						end;
					else
						ur=[]; sr=[]; nr=0;
						if i>2&i<nu-1,
							t_coll(i-1)=collision_time([s(i-2) s(i+1)],...
							[t_start(i-2),t_start(i+1)],...
							[x(i-2),x(i+1)],T+1);
						end;					  
					end;
					s=[s(1:i-2) sr s(i+1:nu-1)];
					hone=ones([1,ns]);
					x=[x(1:i-2) xc*hone x(i+1:nu-1)];
					t_start=[t_start(1:i-2) t*hone t_start(i+1:nu-1)];
					if ns,
						u=[u(:,1:i-1) ur u(:,i+1:nu)];
						t_coll=[t_coll(1:i-1) (T+1)*ones([1,nr]) t_coll(i+1:nu)];
					else
						u=[u(:,1:i-1) u(:,i+2:nu)];
						t_coll=[t_coll(1:i-1) t_coll(i+2:nu)];
					end;
				end;
				[t i]=min(t_coll);
			end;
			n=length(x);
			for i=1:n
				line([x(i),x(i)+s(i)*(T-t_start(i))],[t_start(i) T]);
			end;
			drawnow;
			x=x+s.*(T-t_start);
			hold off;			
		else 
			%% Not saving the graphics
			scn=get(0,'ScreenSize');
			ext=[0.2 0.7];
			pos=[scn(1:2)+ext(1)*(scn(3:4)-scn(1:2)) scn(1:2)+ext(2)*(scn(3:4)-scn(1:2))];
			figure('Position',pos);
			if ~Periodic,
				[u s t_coll,x]=initialize(u0,x0,T,delta);
			else
				[u s t_coll,x]=initialize(u0,x0,T,delta,'periodic',xleft,xright);
			end;
			t_start=zeros(size(x));
			if Periodic,
				xlower=xleft;
				xhigher=xright;
			else
				xlower=min([x(1)+s(1)*T,x(1)]);
				xhigher=max([x(length(x))+s(length(s))*T,x(length(x))]);
				if xlower==xhigher
					xlower=xhigher-0.5;
					xhigher=xlower+1;
				end;
			end;
			axes('XLim',[xlower xhigher],'YLim',[0 T],'Box','on','Drawmode','fast');
			hold on;
			first=1;
			[t i]=min(t_coll);
			ncollision=0;
			while t<T&~isempty(s),
% 				if (~check_states(x,s,t_start,t))
% 					display('not increasing collisions');
% 					ncollision
%				end;
				ncollision=ncollision+1;
				[ncr nu]=size(u);
				nx=length(x);
				if nx~=nu-1,
					[nx nu];
					error('Wrong length of x vs u in FrontTrack');
				end;
				xc=x(i)+s(i)*(t-t_start(i));
% 				if (i>1),
% 					xc2=x(i-1)+(t-t_start(i-1))*s(i-1);
% 					if (abs(xc-xc2)>1e-6)
% 						display(' not same collision point...');
% 					end;
% 				end;
				if first,
					first=0;
					h=line([x(i) xc],[t_start(i) t],'EraseMode','none');
				else
					set(h,'XData',[x(i) xc],'YData',[t_start(i) t]);
				end;
				drawnow;
				if i>1,
					set(h,'XData',[x(i-1) xc],'YData',[t_start(i-1) t]);
				end;
				drawnow;
				if (Periodic&i==2),
					tcoll_left=collision_time([0.0 s(3)], [t t_start(3)], [xleft x(3)],T+1);
					tcoll_right=collision_time([s(nu-2) s(2)], [t_start(nu-2) t],...
						[x(nu-2) xright],T+1);
					x=[xleft x(3:nu-2) xright xright];
					s=[0 s(3:nu-2) s(2) 0];
					t_start=[t t_start(3:nu-2) t t];
					u=[u(:,3) u(:,3:nu-1) u(:,3) u(:,3)];
					t_coll=[T+1 tcoll_left t_coll(4:nu-2) tcoll_right T+1 T+1];
				elseif (Periodic&i==nu-1),
					tcoll_right=collision_time([s(nx-2) 0.0],[t_start(nx-2) t],...
						[x(nx-2) xright],T+1);
					tcoll_left=collision_time([s(nx-1),s(2)],[t t_start(2)],...
						[xleft, x(2)],T+1);
					x=[xleft xleft x(2:nx-2) xright];
					t_start=[t t t_start(2:nx-2) t];
					s=[0 s(nx-1) s(2:nx-2) 0];
					u=[u(:,nu-1) u(:,nu-1) u(:,2:nu-2) u(:,nu-2)];
					t_coll=[T+1,T+1, tcoll_left,t_coll(3:nu-3),tcoll_right,T+1];
				else
					[ur sr]=Riemannsol(u(:,i-1),u(:,i+1),delta);
					[nc nr]=size(ur); 
					ns=length(sr);
					if ns>0,
						if i>2,
							t_coll(i-1)=collision_time([s(i-2) sr(1)],[t_start(i-2),t],...
							[x(i-2),xc],T+1);
							if t_coll(i-1)<t-COLL_TOL,
								[-1 t_coll(i-1) t t-t_coll(i-1)]
								[s(i-2) sr(1) t_start(i-2) t]
								error('not increasing t_coll');
							end;
						end;
						if i<nu-1,
							t_coll(i+1)=collision_time([sr(ns),s(i+1)],[t,t_start(i+1)],...
							[xc,x(i+1)],T+1);
							if t_coll(i+1)<t-COLL_TOL,
								[1 t_coll(i+1) t t-t_coll(i+1)]
								[sr(ns) s(i+1)]
								error('not increasing t_coll');
							end;
						end;
					else
						ur=[]; sr=[]; nr=0;
						if i>2&i<nu-1,
							t_coll(i-1)=collision_time([s(i-2) s(i+1)],...
							[t_start(i-2),t_start(i+1)],...
							[x(i-2),x(i+1)],T+1);
						end;					  
					end;
					s=[s(1:i-2) sr s(i+1:nu-1)];
					hone=ones([1,ns]);
					x=[x(1:i-2) xc*hone x(i+1:nu-1)];
					t_start=[t_start(1:i-2) t*hone t_start(i+1:nu-1)];
					if ns,
						u=[u(:,1:i-1) ur u(:,i+1:nu)];
						t_coll=[t_coll(1:i-1) (T+1)*ones([1,nr]) t_coll(i+1:nu)];
					else
						u=[u(:,1:i-1) u(:,i+2:nu)];
						t_coll=[t_coll(1:i-1) t_coll(i+2:nu)];
					end;
				end;
				[t i]=min(t_coll);
			end;
			n=length(x);
			if first,
				first=0;
				h=line([x(1)+s(1)*(T-t_start(1))],[t_start(1) T],'EraseMode','none');
				drawnow;
			end;	
			for i=1:n
				set(h,'XData',[x(i),x(i)+s(i)*(T-t_start(i))],'YData',[t_start(i) T]);
				drawnow;
			end;		
			x=x+s.*(T-t_start);
			hold off;	
		end;
	elseif polygons,
		%% Using polygons to plot 
		if savegraphics,
			scn=get(0,'ScreenSize');
			ext=[0.2 0.7];
			pos=[scn(1:2)+ext(1)*(scn(3:4)-scn(1:2)) scn(1:2)+ext(2)*(scn(3:4)-scn(1:2))];
			figure('Position',pos);
			% ............For drawing filled polygon ...............
			if cmap_set,
				j=cmap;				
			else
				j=jet(256);				
			end;
			colormap(j); 
			NMAP=length(j);
			% ......................................................
			if ~Periodic,
				[u s t_coll,x]=initialize(u0,x0,T,delta);
			else
				[u s t_coll,x]=initialize(u0,x0,T,delta,'periodic',xleft,xright);
			end;		
		
			cf=colfunc(u);
			pmax=max(cf);pmin=min(cf);
			clear('cf');
			
			t_start=zeros(size(x));
			if Periodic,
				xlower=xleft;
				xhigher=xright;
			else
				xlower=min([x(1)+s(1)*T,x(1)]);
				xhigher=max([x(length(x))+s(length(s))*T,x(length(x))]);
				if xlower==xhigher
					xlower=xhigher-0.5;
					xhigher=xlower+1;
				end;
			end;		
			axes('XLim',[xlower xhigher],'YLim',[0 T],'Box','on','Drawmode','fast');
			hold on;
			[t i]=min(t_coll);
			first=1;
			while t<T&~isempty(s),
				[ncr nu]=size(u);
				nx=length(x);
				if nx~=nu-1,
					[nx nu];
					error('Wrong length of x vs u in FrontTrack');
				end;
				xc=x(i)+s(i)*(t-t_start(i));
				% ............For drawing filled polygon ...............
				xc=x(i)+s(i)*(t-t_start(i));
				xp=[x(i) xc x(i-1)];
				tp=[t_start(i) t t_start(i-1)];
				colind=clindex(u(:,i),pmax,pmin,NMAP);
				patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
				drawnow;
				if i>2,
					xp=[x(i-2) x(i-1) xc];
					tp=[t_start(i-2) t_start(i-1) t];
				else
					xp=[xlower x(i-1) xc xlower];
					tp=[t_start(i-1) t_start(i-1) t t];
				end;
				colind=clindex(u(:,i-1),pmax,pmin,NMAP);
				patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
				if i<nu-2,
					xp=[x(i) x(i+1) xc];
					tp=[t_start(i) t_start(i+1) t];
				else
					xp=[x(i) xhigher xhigher xc];
					tp=[t_start(i) t_start(i) t t];
				end;
				colind=clindex(u(:,i+1),pmax,pmin,NMAP);
				patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
				% ......................................................
				drawnow;
				if (Periodic&i==2),
					xp=[x(nx-1) x(nx) x(nx)];
					tp=[t_start(nx-1) t_start(nx) t];
					colind=clindex(u(:,nx),pmax,pmin,NMAP);
					patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
                    drawnow;			
					tcoll_left=collision_time([0.0 s(3)], [t t_start(3)], [xleft x(3)],T+1);
					tcoll_right=collision_time([s(nu-2) s(2)], [t_start(nu-2) t],...
						[x(nu-2) xright],T+1);
					x=[xleft x(3:nu-2) xright xright];
					s=[0 s(3:nu-2) s(2) 0];
					t_start=[t t_start(3:nu-2) t t];
					u=[u(:,3) u(:,3:nu-1) u(:,3) u(:,3)];
					t_coll=[T tcoll_left t_coll(4:nu-2) tcoll_right T+1 T+1];
				elseif (Periodic&i==nu-1),
					xp=[x(1) x(1) x(2)];
					tp=[t_start(1) t t_start(2)];
					colind=clindex(u(:,2),pmax,pmin,NMAP);
					patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
					tcoll_right=collision_time([s(nx-2) 0.0],[t_start(nx-2) t],...
						[x(nx-2) xright],T+1);
					tcoll_left=collision_time([s(nx-1),s(2)],[t t_start(2)],...
						[xleft, x(2)],T+1);
					x=[xleft xleft x(2:nx-2) xright];
					t_start=[t t t_start(2:nx-2) t];
					s=[0 s(nx-1) s(2:nx-2) 0];
					u=[u(:,nu-1) u(:,nu-1) u(:,2:nu-2) u(:,nu-2)];
					t_coll=[T+1,T+1, tcoll_left,t_coll(3:nu-3),tcoll_right,T+1];
				else
					[ur sr]=Riemannsol(u(:,i-1),u(:,i+1),delta);
					[nc nr]=size(ur); 
					ns=length(sr);
					if ns>0,
						if i>2,
							t_coll(i-1)=collision_time([s(i-2) sr(1)],[t_start(i-2),t],...
							[x(i-2),xc],T+1);
							if t_coll(i-1)<t-COLL_TOL,
								[-1 t_coll(i-1) t t-t_coll(i-1)]
								[s(i-2) sr(1) t_start(i-2) t]
								error('not increasing t_coll');
							end;
						end;
						if i<nu-1,
							t_coll(i+1)=collision_time([sr(ns),s(i+1)],[t,t_start(i+1)],...
							[xc,x(i+1)],T+1);
							if t_coll(i+1)<t-COLL_TOL,
								[1 t_coll(i+1) t t-t_coll(i+1)]
								[sr(ns) s(i+1)]
								error('not increasing t_coll');
							end;
						end;
					else
						ur=[]; sr=[]; nr=0;
						if i>2&i<nu-1,
							t_coll(i-1)=collision_time([s(i-2) s(i+1)],...
							[t_start(i-2),t_start(i+1)],...
							[x(i-2),x(i+1)],T+1);
						end;					  
					end;
					s=[s(1:i-2) sr s(i+1:nu-1)];
					hone=ones([1,ns]);
					x=[x(1:i-2) xc*hone x(i+1:nu-1)];
					t_start=[t_start(1:i-2) t*hone t_start(i+1:nu-1)];
					if ns,
						u=[u(:,1:i-1) ur u(:,i+1:nu)];
						t_coll=[t_coll(1:i-1) (T+1)*ones([1,nr]) t_coll(i+1:nu)];
					else
						u=[u(:,1:i-1) u(:,i+2:nu)];
						t_coll=[t_coll(1:i-1) t_coll(i+2:nu)];
					end;
				end;
				[t i]=min(t_coll);
			end;
			% ...............For drawing filled polygon ..............
			n=length(x);
			xp=[xlower xlower];
			tp=[t_start(1) T]; 
			for i=1:n
				xp2=[x(i)+s(i)*(T-t_start(i)) x(i)];
				tp2=[T t_start(i)];
				colind=clindex(u(:,i),pmax,pmin,NMAP);
				patch([xp xp2],[tp tp2],j(colind,:),'EdgeColor','none','EraseMode','none');
				xp=xp2(2:-1:1); tp=tp2(2:-1:1);
			end;
			drawnow;
			i=n;
			xp2=[xhigher xhigher]; tp2=[T t_start(i)];
			colind=clindex(u(:,i+1),pmax,pmin,NMAP);
			patch([xp xp2],[tp tp2],j(colind,:),'EdgeColor','none','EraseMode','none');
			% ..........................................................
			x=x+s.*(T-t_start); 
			hold off;
		else
			%% Polygons without graphics save ...
			scn=get(0,'ScreenSize');
			ext=[0.2 0.7];
			pos=[scn(1:2)+ext(1)*(scn(3:4)-scn(1:2)) scn(1:2)+ext(2)*(scn(3:4)-scn(1:2))];
			figure('Position',pos);
			% ............For drawing filled polygon ...............
			if cmap_set,
				j=cmap;				
			else
				j=jet(256);				
			end;
			colormap(j); 
			NMAP=length(j);
			% ......................................................
			if ~Periodic,
				[u s t_coll,x]=initialize(u0,x0,T,delta);
			else
				[u s t_coll,x]=initialize(u0,x0,T,delta,'periodic',xleft,xright);
			end;		
		
			cf=colfunc(u);
			pmax=max(cf);pmin=min(cf);
			clear('cf');
			
			t_start=zeros(size(x));
			if Periodic,
				xlower=xleft;
				xhigher=xright;
			else
				xlower=min([x(1)+s(1)*T,x(1)]);
				xhigher=max([x(length(x))+s(length(s))*T,x(length(x))]);
				if xlower==xhigher
					xlower=xhigher-0.5;
					xhigher=xlower+1;
				end;
			end;		
			axes('XLim',[xlower xhigher],'YLim',[0 T],'Box','on','Drawmode','fast');
			hold on;
			[t i]=min(t_coll);
			first=1;
			while t<T&~isempty(s),
				[ncr nu]=size(u);
				nx=length(x);
				if nx~=nu-1,
					[nx nu];
					error('Wrong length of x vs u in FrontTrack');
				end;
				xc=x(i)+s(i)*(t-t_start(i));
				% ............For drawing filled polygon ...............
				xc=x(i)+s(i)*(t-t_start(i));
				xp=[x(i) xc x(i-1)];
				tp=[t_start(i) t t_start(i-1)];
				colind=clindex(u(:,i),pmax,pmin,NMAP);
				if first,
					h=patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
					first=0;
				else
					set(h,'XData',xp,'YData',tp,'FaceColor',j(colind,:));
				end;
				drawnow;
				if i>2,
					xp=[x(i-2) x(i-1) xc];
					tp=[t_start(i-2) t_start(i-1) t];
				else
					xp=[xlower x(i-1) xc xlower];
					tp=[t_start(i-1) t_start(i-1) t t];
				end;
				colind=clindex(u(:,i-1),pmax,pmin,NMAP);
				set(h,'XData',xp,'YData',tp,'FaceColor',j(colind,:));
				drawnow;
				if i<nu-2,
					xp=[x(i) x(i+1) xc];
					tp=[t_start(i) t_start(i+1) t];
				else
					xp=[x(i) xhigher xhigher xc];
					tp=[t_start(i) t_start(i) t t];
				end;
				colind=clindex(u(:,i+1),pmax,pmin,NMAP);
				set(h,'XData',xp,'YData',tp,'FaceColor',j(colind,:));
				drawnow;
				if (Periodic&i==2),
					xp=[x(nx-1) x(nx) x(nx)];
					tp=[t_start(nx-1) t_start(nx) t];
					colind=clindex(u(:,nx),pmax,pmin,NMAP);
					if first,
						first=0,
						h=patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
					else
						set(h,'XData',xp,'YData',tp,'FaceColor',j(colind,:));
					end;				
					tcoll_left=collision_time([0.0 s(3)], [t t_start(3)], [xleft x(3)],T+1);
					tcoll_right=collision_time([s(nu-2) s(2)], [t_start(nu-2) t],...
						[x(nu-2) xright],T+1);
					x=[xleft x(3:nu-2) xright xright];
					s=[0 s(3:nu-2) s(2) 0];
					t_start=[t t_start(3:nu-2) t t];
					u=[u(:,3) u(:,3:nu-1) u(:,3) u(:,3)];
					t_coll=[T tcoll_left t_coll(4:nu-2) tcoll_right T T];
				elseif (Periodic&i==nu-1),
					xp=[x(1) x(1) x(2)];
					tp=[t_start(1) t t_start(2)];
					colind=clindex(u(:,2),pmax,pmin,NMAP);
					if first,
						first=0,
						h=patch(xp,tp,j(colind,:),'EdgeColor','none','EraseMode','none');
					else
						set(h,'XData',xp,'YData',tp,'FaceColor',j(colind,:));
					end;
					tcoll_right=collision_time([s(nx-2) 0.0],[t_start(nx-2) t],...
						[x(nx-2) xright],T+1);
					tcoll_left=collision_time([s(nx-1),s(2)],[t t_start(2)],...
					[xleft, x(2)],T+1);
					x=[xleft xleft x(2:nx-2) xright];
					t_start=[t t t_start(2:nx-2) t];
					s=[0 s(nx-1) s(2:nx-2) 0];
					u=[u(:,nu-1) u(:,nu-1) u(:,2:nu-2) u(:,nu-2)];
					t_coll=[T,T, tcoll_left,t_coll(3:nu-3),tcoll_right,T];
				else
					[ur sr]=Riemannsol(u(:,i-1),u(:,i+1),delta);
					[nc nr]=size(ur); 
					ns=length(sr);
					if ns>0,
						if i>2,
							t_coll(i-1)=collision_time([s(i-2) sr(1)],[t_start(i-2),t],...
							[x(i-2),xc],T+1);
							if t_coll(i-1)<t-COLL_TOL,
								[-1 t_coll(i-1) t t-t_coll(i-1)]
								[s(i-2) sr(1) t_start(i-2) t]
								error('not increasing t_coll');
							end;
						end;
						if i<nu-1,
							t_coll(i+1)=collision_time([sr(ns),s(i+1)],[t,t_start(i+1)],...
							[xc,x(i+1)],T+1);
							if t_coll(i+1)<t-COLL_TOL,
								[1 t_coll(i+1) t t-t_coll(i+1)]
								[sr(ns) s(i+1)]
								error('not increasing t_coll');
							end;
						end;
					else
						ur=[]; sr=[]; nr=0;
						if i>2&i<nu-1,
							t_coll(i-1)=collision_time([s(i-2) s(i+1)],...
							[t_start(i-2),t_start(i+1)],...
							[x(i-2),x(i+1)],T+1);
						end;					  
					end;
					s=[s(1:i-2) sr s(i+1:nu-1)];
					hone=ones([1,ns]);
					x=[x(1:i-2) xc*hone x(i+1:nu-1)];
					t_start=[t_start(1:i-2) t*hone t_start(i+1:nu-1)];
					if ns,
						u=[u(:,1:i-1) ur u(:,i+1:nu)];
						t_coll=[t_coll(1:i-1) (T+1)*ones([1,nr]) t_coll(i+1:nu)];
					else
						u=[u(:,1:i-1) u(:,i+2:nu)];
						t_coll=[t_coll(1:i-1) t_coll(i+2:nu)];
					end;
				end;
				[t i]=min(t_coll);
			end;
			% ...............For drawing filled polygon ..............
			n=length(x);
			xp=[xlower xlower];
			tp=[t_start(1) T]; 
			for i=1:n
				xp2=[x(i)+s(i)*(T-t_start(i)) x(i)];
				tp2=[T t_start(i)];
				colind=clindex(u(:,i),pmax,pmin,NMAP);
				if first,
					first=0;
					h=patch([xp xp2],[tp tp2],'Facecolor',j(colind,:),'EraseMode','none');
				else
					set(h,'Xdata',[xp xp2],'YData',[tp tp2],'FaceColor',j(colind,:));
				end;
				xp=xp2(2:-1:1); tp=tp2(2:-1:1);
			end;
			drawnow;
			i=n;
			xp2=[xhigher xhigher]; tp2=[T t_start(i)];
			colind=clindex(u(:,i+1),pmax,pmin,NMAP);
			set(h,'Xdata',[xp xp2],'YData',[tp tp2],'FaceColor',j(colind,:));
			drawnow;
			x=x+s.*(T-t_start); 
			hold off;
		end;
	end;
end;
%%
