function [u,s,t_coll,x]=initialize(u0,x0,T,delta,varargin)
%
%	Initializes the fronts to use by fronttracking. The initial values are 
% 	in the vector u0(1:#components,1:n) and the positions of the 
%   discontinuities in x0(1:n-1). delta gives the accuracy of the Riemann solver
%	and T is the time to front track, this is needed to set the collision time.
%
Periodic=0;
l=length(varargin);
if l>0,
	if strcmp(varargin(1),'periodic')
		Periodic=1;
		bndr=varargin(2);
		xleft=bndr{1};
		bndr=varargin(3);
		xright=bndr{1};
	else
		error('                  Unknown option in initialize');
	end;
end;
[nc n]=size(u0);
i=1;
while (max(abs(u0(:,i)-u0(:,i+1)))<delta*0.001)&(i<n),
	i=i+1;
end;
ni=i;
[ur s]=Riemannsol(u0(:,ni),u0(:,ni+1),delta);
u=[u0(:,ni) ur u0(:,ni+1)];
[nc nlast]=size(u);
x=x0(ni)*ones(size(s));
t_coll=(T+1)*ones([1,nlast]);
for i=ni+1:n-1
	[ur sri]=Riemannsol(u(:,nlast),u0(:,i+1),delta);
	nx=length(x);
	if ~isempty(sri),
		u=[u ur u0(:,i+1)];
		s=[s sri];
		t_coll(nlast)=collision_time([s(nlast-1) sri(1)],[0,0],[x(nx),x0(i)]);
		ns=length(sri);
		t_coll=[t_coll (T+1)*ones([1,ns])];
		x=[x x0(i)*ones([1,ns])];
	end;
	[nc nlast]=size(u);
end;
if Periodic,
	if ~all(u(:,1)==u(:,nlast)),
		error('                  Not correct periodic initial data');
	end;
	if x(1)<=xleft,
		error('                  Not proper left boundary for periodic');
	end;
	if x(nlast-1)>=xright,
		error('                  Not proper right boundary for periodic');
	end;
	tcoll_left=collision_time([0.0 s(1)],[0 0],[xleft x(1)]);
	tcoll_right=collision_time([s(nlast-1) 0.0], [0 0], [x(nlast-1) xright]);
	x=[xleft x xright];
	s=[0 s 0];
	t_coll=[T+1 tcoll_left t_coll(2:nlast-1) tcoll_right T+1];
	u=[u(:,1) u u(:,nlast)];
end;


