clear
clc
close all

N_1 = 16;
N_2 = 16;
dt = 0.1;

N_car = N_1 + N_2;

index_car = randperm(N_car);

index_1 = index_car(1: N_1) ;
index_2 = [index_car(N_1+1:N_1+N_2) ];

N_T = 2001;
car_length = 5;

desired_space = 25; 
delta_u = desired_space + car_length;

v_avg = 25;

k_d = zeros(1, N_car);
k_v = k_d;

k_d(index_1)= 0.1 + 0.1*(rand(size(index_1)) - 0.5);
k_v(index_1)= 0.3 + 0.2*(rand(size(index_1)) - 0.5);

k_d(index_2)= 0.3 + 0.2*(rand(size(index_2)) - 0.5);
k_v(index_2)= 0.1 + 0.1*(rand(size(index_2)) - 0.5);

u0 = 20*rand(size(index_car)) - 10;
v0 = 10*rand(size(index_car)) - 5;
d2u = zeros(size(u0));
d2v = d2u;
du = d2u;
dv = d2v;

U=zeros(N_T, N_car );
V = U;

U(1, :) = u0;
V(1,:) = v0;

U_shift = (2 : N_car+2)*delta_u;

u =  U_shift + [0, u0];

Time = [0:N_T-1]*dt;

Y = zeros(1, N_car);

a_min = -5;
a_max = 5;
v_max = v_avg;
v_min = -19.444;

Total_length = (N_car + 4)*delta_u;
dx_hist = Total_length/N_car;


figure(1);
set(gcf, 'position',[60,150,1300,150])
subplot('position',[0.03, 0.3, 0.94, 0.4])
plot(u(index_1+1), Y(index_1),'sr', 'MarkerSize',14,'MarkerFaceColor','w','lineWidth',3)
hold on
plot(u(index_2+1), Y(index_2),'sk', 'MarkerSize',14,'MarkerFaceColor','w', 'lineWidth',3)
plot(u(1), Y(1),'sk', 'MarkerSize',14,'MarkerFaceColor','k', 'lineWidth',3)
title('Inhomogeneous BCC traffic flow (time: 0.00 sec.)')

xlim([0, Total_length]);
ylim([-0.5, 0.5]);
set(gca,'ytick',[], 'FontSize', 16)
set(gca,'xtick', [0:4:N_car]*dx_hist, 'FontSize',20);
hold off

M(1) = getframe(gcf);

U0 = zeros(N_T,1);
V0 = U0;

for i = 2:N_T

    if i>1000 && i<=1020
        a=5;
        v0 = V0(i-1) + a*dt;
    elseif i>1020 && i<=1040
        a = -5;
        v0 = V0(i-1) + a*dt;
    else
        v0 = 0;
        a = 0;
    end
     u0 = U0(i-1) + v0*dt + a*dt^2/2;
     U0(i) = u0;
    V0(i) = v0;
   
    d2u(1) =  u0 + U(i-1,2) - 2*U(i-1,1);
    d2u(end) = U(i-1,end-1) - U(i-1,end);
    d2u(2:end-1) = U(i-1,1:end-2) + U(i-1,3:end) - 2*U(i-1,2:end-1);
    d2u =  k_d .* d2u;
    
    d2v(1) =  v0 + V(i-1,2) - 2*V(i-1,1);
    d2v(end) = V(i-1,end-1) - V(i-1,end);
    d2v(2:end-1) = V(i-1,1:end-2) + V(i-1,3:end) - 2*V(i-1,2:end-1);
    d2v =  k_v .* d2v;
    
    a = max(a_min,min(a_max,d2u+d2v));
    v = max(v_min,min(v_max,V(i-1,:) + a*dt)) ;
    u = U(i-1,:) + v*dt + a*dt^2/2;
    
    

    V(i,:) = v;
    U(i,:) = u;
    
    U_shift = mod(U_shift - v_avg * dt, Total_length);
    
    u =  U_shift + [U0(i), U(i,:)];

figure(1);
set(gcf, 'position',[60,150,1300,150])
subplot('position',[0.03, 0.3, 0.94, 0.4])
plot(u(index_1+1), Y(index_1),'sr', 'MarkerSize',14,'MarkerFaceColor','w','lineWidth',3)
hold on
plot(u(index_2+1), Y(index_2),'sk', 'MarkerSize',14,'MarkerFaceColor','w', 'lineWidth',3)
plot(u(1), Y(1),'sk', 'MarkerSize',14,'MarkerFaceColor','k', 'lineWidth',3)
title(sprintf('Inhomogeneous BCC traffic flow (time: %3.1f sec.)',(i-1)*dt));

xlim([0, Total_length]);
ylim([-0.5, 0.5]);
set(gca,'ytick',[], 'FontSize', 16)
set(gca,'xtick', [0:4:N_car]*dx_hist, 'FontSize',20);
hold off

M(i) = getframe(gcf);
   
end


