function U = wvtraffn(x,t,v_max,u_max,u_0,L)
%WVTRAFFN  u = wvtraffn(x,t,v_max,u_max,u_0,L) returns the value u(x,t) of the
%          solution of the traffic flow model 
%          
%          u_t + v_max(1 - 2u/u_max) u_x = 0
%          u(x,0) = u_0    x < -L
%                   u_max  -L < x < 0
%                   0      x > 0
%
%           

%Created: 5/95 by R. Knobel


t_0 = u_max*L / (v_max*(u_max-u_0));

if t==0,

  U = u_0*(-x-L >= 0) + u_max*(x+L-eps >= 0).*(-x >= 0);

elseif (t > 0 & t <= t_0),

  x0 = -v_max*u_0*t/u_max - L;
  x1 = -v_max*t;
  x2 =  v_max*t;
  U  = u_0*(x0-x >= 0) + u_max*(x-x0-eps >= 0).*(x1-x >= 0) +...
      (u_max/2)*(1- (1/(v_max*t))*x).*(x-x1-eps >= 0).*(x2-x >= 0);

elseif (t > t_0),

  x0 = v_max*(1 - 2*u_0/u_max)*t - 2*sqrt( L*v_max*t*(1-u_0/u_max) );
  x2 = v_max*t;
  U  = u_0*(x0-x >= 0) + ...
       (u_max/2)*(1- (1/(v_max*t))*x).*(x-x0-eps >= 0).*(x2-x >= 0);

end;