%% This code mainly comes from https://blog.csdn.net/weixin_43102634/article/details/102996254?ops_request_misc=&request_id=&biz_id=102&utm_term=%E4%B8%80%E7%BB%B4%E5%85%83%E8%83%9E%E8%87%AA%E5%8A%A8%E6%9C%BAmatlab&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduweb~default-9-102996254.first_rank_v2_pc_rank_v29&spm=1018.2226.3001.4187
% Among them, lines 37-90 realize the modeling of cellular automata
% Lines 44-74 add the producer-consumer model
% @Shao Yi, 2022/04/04


clc
clear;
%%
%build the GUI
%define the plot button
plotbutton=uicontrol('style','pushbutton','string','Run', 'fontsize',12,'position',[100,400,50,20],'callback','run=1;');
%define the stop button
erasebutton=uicontrol('style','pushbutton','string','Stop','fontsize',12,'position',[200,400,50,20],'callback','freeze=1;');
%define the Quit button
quitbutton=uicontrol('style','pushbutton','string','Quit','fontsize',12,'position',[300,400,50,20],'callback','stop=1;close;');
number=uicontrol('style','text','string','1','fontsize',12,'position',[20,400,50,20]);

%%
%CAsetup
n=1000; %Data initialization
z=zeros(1,n); %Cell number
z=roadstart(z,200); %The road state was initialized, and 200 vehicles were randomly distributed on the road section
cells=z;
vmax=5; %Maximum speed
v=speedstart(cells,vmax); 
x=1;
memor_cells=zeros(1000,n);
memor_v=zeros(1000,n);
imh=imshow(cells); %Initialize image white with car, black with empty cell
set(imh,'erasemode','none')
axis equal
axis tight
stop=0; %wait for a quit button push
run=0; %wait for a draw
freeze=0; %wait for a freeze 
while (stop==0 && x<1000)
if(run==1)
      %Boundary condition processing, search first and last vehicles, control entry and exit, use opening conditions
      a=searchleadcar(cells);
      b=searchlastcar(cells);
%       [cells,v]=border_control(cells,a,b,v,vmax);
      i=searchleadcar(cells); %Search the location of the first train
      index=randperm(i);
      for j=1:index
      if (i-j+1==n)
      [z,v]=leadcarupdate(z,v);
            continue;
          else
      
      if cells(i-j+1)==0 
           continue;
      else
          v(i-j+1)=min(v(i-j+1)+1,vmax); %accelerate
  
      k=searchfrontcar((i-j+1),cells); 
      if(k==0) 
      d=n-(i-j+1);
      else
          d=k-(i-j+1)-1;
      end
      v(i-j+1)=min(v(i-j+1),d);%decelerate
      %Random slowing down
      v(i-j+1)=randslow(v(i-j+1));
      new_v=v(i-j+1);
      %Update vehicle Location
      z(i-j+1)=0;
      z(i-j+1+new_v)=1;
      %Update speed
      v(i-j+1)=0;
      v(i-j+1+new_v)=new_v;
      end
      end
      end
      cells=z;
      memor_cells(x,:)=cells; %Record speed and vehicle position
      memor_v(x,:)=v;
      x=x+1;
      set(imh,'cdata',cells) %Update the image
      %update the step number diaplay
pause(0.0001);
stepnumber=1+str2num(get(number,'string'));
set(number,'string',num2str(stepnumber)) 
end
if (freeze==1)
run=0;
freeze=0;
end
drawnow 
end
%% Figure 2
figure(2)
for l=1:1:200
for k=500:1:1000
if memor_cells(l,k)>0
    plot(k,l,'k.');
hold on;
end
end
end
xlabel('spatial location');
ylabel('time(s)');
title('Space-time figure');
  % Figure 3
      for i=1:1:1000
      density(i)=sum(memor_cells(i,:)>0)/1000;
      flow(i)=sum(memor_v(i,:))/1000;
      end
      figure(3)
      plot(density,flow,'k.');
      title('Flow-density diagram')
      xlabel('density')
      ylabel('flow')
 %%
  %Figure 4
    figure(4)
       imagesc(memor_v);
       colorbar;
    hold on;
    xlabel('spatial location');
    ylabel('time(s)');
    title('Thermodynamic diagram'); 
 %%
 
 % Function: speedstart.m program code     
function [v_matixcells]=speedstart(matrix_cells,vmax)
    %Road initial state Vehicle speed initializes
    v_matixcells=zeros(1,length(matrix_cells));
    for i=1:length(matrix_cells)
        if matrix_cells(i)~=0
        v_matixcells(i)=round(vmax* rand(1));
        end
    end
end

%Function: searchleadcar.m program code
function [location_leadcar]=searchleadcar(matrix_cells)
    i=length(matrix_cells);
    for j=1:i
        if matrix_cells(i-j+1)~=0
            location_leadcar=i-j+1;
            break;
        else
            location_leadcar=0;
        end
    end
end

%Function: leadcarupdate.m program code
function [new_matrix_cells,new_v]=leadcarupdate(matrix_cells,v)

    n=length(matrix_cells);
    if v(n)~=0
        matrix_cells(n)=0;
        v(n)=0;
    end
    new_matrix_cells=matrix_cells;
    new_v=v;
end

%Function: randslow.m program code
function[new_v]=randslow(v)
    p=0.3; %The probability of slowing down
    rand('state',sum(100*clock)*rand(1));
    p_rand=rand; %Generate random probability
    if p_rand<=p
        v=max(v-1,0);
    end
    new_v=v;
end
 

%Function: roadstart.m program code
function [matrix_cells_start]=roadstart(matrix_cells,n)
%The initializing state of vehicles on the road, cell matrix is randomly 0 or 1, matrix_cells initial matrix, n initial number of vehicles
    k=length(matrix_cells);
    z=round(k*rand(1,n));
    for i=1:n
        j=z(i);
        if j==0
            matrix_cells(j)=0;
        else
            matrix_cells(j)=1;
        end
    end
    matrix_cells_start=matrix_cells;
end

% Function：searchfrontcar.m program code
function [location_frontcar]=searchfrontcar(current_location,matrix_cells)
    i=length(matrix_cells);
    if current_location==i
        location_frontcar=0;
    else
        for j=current_location+1:i
            if matrix_cells(j)~=0
            location_frontcar=j;
            break;
            else
            location_frontcar=0;
            end
         end
    end
end


%Function: searchlastcar.m program code
function [location_lastcar]=searchlastcar(matrix_cells)
%Search for last car location
for i=1:length(matrix_cells)
    if matrix_cells(i)~=0
        location_lastcar=i;
        break;
      else 
          location_lastcar=length(matrix_cells );
    end
end
end