clear all
close all
clc

set(gca,'FontSize',18)
linewidth = 2;

x = -5:0.01:5;
rho = 1/2*(1+tanh(x));
c = 1 - 2*rho;

set(gca,'FontSize',18)
plot(x,rho,'Linewidth',linewidth)
xlabel('x')
ylabel('\rho')
print -djpeg traffic_example2_rho


figure
set(gca,'FontSize',18)
plot(x,c,'Linewidth',linewidth)
xlabel('x')
ylabel('c')
print -djpeg traffic_example2_c

% first draw the characteristics

x_0 = [-10:-1,-0.5,-0.2,-0.1];x_0 = [x_0,0,-x_0];
rho_0 = 1/2*(1+tanh(x_0));
c_0 = 1 - 2*rho_0;
t = 0:0.01:5;

figure

for j = 1:length(x_0)
    set(gca,'FontSize',18)
    plot(x_0(j) + c_0(j)*t,t,'Linewidth',linewidth)
    axis([-5,5,0,5])
    hold on
end
hold off
xlabel('x')
ylabel('t')
print -djpeg traffic_example2_characteristics

figure

x = -10:0.01:10;
rho = 1/2*(1+1*tanh(x));
c = 1 - 2*rho;

t_range=0:0.2:5;
delete traffic_example2.gif    %have to delete old one so new one shows up
for j=1:length(t_range)
    set(gca,'FontSize',18)
    t=t_range(j);
    plot(x + c*t,rho,'r--','LineWidth',linewidth)
    hold on
    % now find the positive part
    [y,i] = find((x + c*t)>0);
    plot(x(1:i(1)) + c(1:i(1))*t,rho(1:i(1)),'LineWidth',linewidth)
    [y,i2] = find((x + c*t)<0);
    plot(x(i2(end):end) + c(i2(end):end)*t,rho(i2(end):end),'LineWidth',linewidth)
    plot([x(i(1)),x(i2(end))] + [c(i(1)),c(i2(end))]*t,[rho(i(1)),rho(i2(end))],'LineWidth',linewidth)
    grid on
    hold off
    axis([-5,5,0,1.1])
    title(['t=',num2str(t_range(j) - 0.01,'%1.1f')])
    ylabel('\rho(x,t)')
    xlabel('x')
    pause(0.1)
    gif_add_frame(gcf,'traffic_example2.gif',1);
end