
% Solution of the traffic flow problem
close all;
clear all;

disp('  ')
disp('Enter 1 for decreasing profile  ')
disp('Enter 2 for increasing profile  ')

m = input( 'enter the choice of initial data:  ')

if m==1 % g(x)
    
    %ymiddle=beta*(1.0-x.^2.*(3.0-2.0*x));
    f_initial = @(x,beta) beta*(x<0)+( (x<1.001)-(x<0) ).*beta.*(1.0-x.^2.*(3.0-2.0*x));
    
elseif m==2 % f(x)
    
    f_initial=@(x,beta) beta*(x>1)+( (x<1.001)-(x<0) ).*beta.*x.^2.*(3.0-2.0*x);
end;
%beta=1.0;
%ymiddle=beta*x.^2.*(3.0-2.0*x);
%y=beta*(x>1)+( (x<1.001)-(x<0) ).*ymiddle;
%% Characteristic speed c(u)=beta-2u

char_speed=@(u,beta)beta-2*u;
%% Characteristic curves starting from x0: x(t)=x0+c(f(x0))*t

char_lines=@(x0,t,beta) x0+char_speed(f_initial(x0,beta),beta).*t;
    
x0=[-1:.1:2];
t=[0:.3:1];
beta=1;

for i=1:length(x0)
    plot(char_lines(x0(i),t,beta),t) 
    hold on
end;
hold off

xlabel('x axis','FontSize',14);
ylabel('t axis','FontSize',14);
title('Characteristics')   
%% Plot the solution at some times t=t1, t2, t3, t4

t = input(' Enter the times in the form  [t1 t2 t3 t4 ]  ')

x = -1:.05:2;
y=f_initial(x,beta);

x1 = x + t(1)*char_speed(y,beta);
x2 = x + t(2)*char_speed(y,beta);
x3 = x + t(3)*char_speed(y,beta);
x4 = x + t(4)*char_speed(y,beta);

figure

plot(x,f_initial(x,beta), 'r')
   
hold on
plot(x1,y,'b-o')   
plot(x2,y,'y-+')
plot(x3,y,'k-.')
plot(x4,y,'g--') 

hold off

say1=['t=' num2str(t(1))];
say2=['t=' num2str(t(2))];
say3=['t=' num2str(t(3))];
say4=['t=' num2str(t(4))];

legend('t=0',say1,say2,say3,say4,'Location', 'Southeast' )

%% Plot the solution surface

%Plot the solution in 3D u=f(x-c(f(x)).*t);
u_sol=@(x,t,beta)f_initial(x-char_speed(f_initial(x,beta),beta).*t,beta);
xx=-1:.05:2;
t=0:.05:1;        % t-axis
[X,Y]=meshgrid(xx,t);      
Z=u_sol(X,Y,beta);

figure
surf(X,Y,Z,'LineStyle','none');
colorbar
xlabel('x','FontSize',14);
ylabel('t','FontSize',14);
zlabel('u','FontSize',14);

figure
pcolor(X,Y,Z);