%  Traffic Simulation using the Nagel-Schreckenberg Model
%  with optional VDR(Velocity Dependent Randomization)
%============================================================%
%  Consider a circlular, one lane road 'road_length' units   %
%  long.  The road is populated with 'num_vehicles' cars.    %
%  Cars are not allowed to collide with or pass each other.  %
%  Cars start with an initial velocity and run through a     %
%  four step simulation:  1)  Cars accelerate by one unit,   %
%  2)  Collision is checked for and avoided, 3)  The cars    %
%  are randomly slowed down(this is the Monte Carlo step),   %
%  4)  Each car advances to it's new position (time step).   %
%============================================================%
%  ABW 2015  %
%============%
clc;

vdr_switch = 'off';  %  Turn velocity dependent randomization on/off

%  Initialize model parameters
road_length     = 200;          %  Length of road
iter_max        = 400;          %  Number of time steps
num_vehicles    = 50;           %  Number of vehicles
max_speed       = 5;            %  Speed limit
p_slow_mov      = 0.01;         %  Probability of a moving car slowing
p_slow_std      = 0.5;          %  Probability of a standing car slowing

%  Initialize vectors
d = zeros(1,num_vehicles);      %  Holds distance between cars
vel = zeros(1,num_vehicles);    %  Holds velocity
pslow = zeros(1,num_vehicles);  %  Holds slowing probability

%  Initialize posn matrix
%  0 represents empty road, 1 represents road occupied by one car
posn = zeros(iter_max,road_length);

%  Uniformly distribute num_vehicles on road
icount = 0;
for i = 1:road_length
    if mod(i,floor(road_length/num_vehicles)) == 0 && icount < num_vehicles
        posn(1,i) = 1;
        icount = icount + 1;
    end
end

%  Set starting velocity to random number or set to zero
for i = 1:num_vehicles
   %vel(i) = randi(max_speed+1)-1; %  Random number in range [0,max_speed]
   vel(i) = 0;
end

for tstep = 1:iter_max-1
    
    %  Set up dynamic randomization for step 3 (VDR)
    for i = 1:num_vehicles
        if vel(i) > 0
           pslow(i) = p_slow_mov;
        elseif vel(i) < 1
           pslow(i) = p_slow_std;
        end
    end
    
    %  Begin Step 1, Acceleration
    for i = 1:num_vehicles
        if vel(i) < max_speed
            vel(i) = vel(i) + 1;
        end
    end
    
    %  Begin Step 2, Braking
    icount = 0;
    dcount = 0;
    for i = 1:road_length
        if i == 1 && posn(tstep,i) == 1
            icount = 1;
            i_current = 1;
        end
        if i > 1 && posn(tstep,i) == 1 && icount < num_vehicles
            if icount > 0
                d(icount) = dcount;  % Write distance for current car
                dcount = 0;  %  Reset distance measurement
                icount = icount + 1;  %  Move on to next car
            elseif icount < 1
                dcount = 0;
                i_current = i;
                icount = icount + 1;
            end
        end
        if posn(tstep,i) == 0 && icount < num_vehicles
            dcount = dcount + 1;  %  On empty road, keep counting
        end
        if icount == num_vehicles
            d(num_vehicles) = (road_length - i) + (i_current - 1);
            icount = 0;
        end
    end
    for i = 1:num_vehicles
        if vel(i) > d(i)
            vel(i) = d(i);
        end
    end
    
    %  Begin Step 3, Random Slowing
    %  If VDR is turned off, pslow = p_slow_std
    if strcmp(vdr_switch,'on') == 1
        for i = 1:num_vehicles
            if rand < pslow(i) && vel(i) > 0
                vel(i) = vel(i) - 1; 
            end
        end
    else
        for i = 1:num_vehicles
            if rand < p_slow_std && vel(i) > 0
              vel(i) = vel(i) - 1; 
            end
        end
    end
    
    %  Begin Step 4, Motion
    icount = 0;
    new_posn = 0;
    for i = 1:road_length
        if posn(tstep,i) == 1
            new_posn = i + vel(icount + 1);
            if new_posn > road_length
                new_posn = new_posn - road_length;
            end
            posn(tstep+1,new_posn) = 1;
            icount = icount + 1;
        end
    end
end

%  Sparse Matrix Plot with spy command
spy(posn,'r.',1);
xlabel('<-- Distance -->')
ylabel('<-- Time Step -->')
title('Sparse Matrix Plot')