clear all;clc
density(:,1)=[30,30,30,30,30,30,30,30,30,30,30,30,30,30]';
velocity(:,1)=[50,50,50,50,50,50,50,50,50,50,50,50,50,50]';
flow(:,1)=[1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,]';
velocity_free=80;
density_jam=80;
density_d=36.75;
ramp_min=200;
T=0.00417;     
L=0.5;
tao=0.01;                                                    
velocity_free=80;                                                   
gamma=35;                                                                                                          
k_constant=13;                                                
alpha=0.95;                                                   
density_jam=80;
fb_gain=30;                                                                %queue_max=;%queue_min=; To be tried later.
ramp_max=1800;
ramp_min=340;

for i=1:14
    for k=1:600
        off_ramp(i,k)=0;
        on_ramp(i,k)=0;
    end
end
On_ramp3=zeros(1,600);
On_ramp10=zeros(1,600);
Demand3=zeros(1,600);
Demand10=zeros(1,600);
Off_ramp=zeros(1,600);
Flow=1500;                                                                 %Flow=zeros(1,600);
queue3=zeros(1,601);
queue10=zeros(1,601);
on_ramp3=zeros(1,600);
on_ramp10=zeros(1,600);

for i=1:600
    if i<200
        Demand3(i)=100+0.5*i;
        Demand10(i)=160+0.5*i;
        Off_ramp(i)=150+0.9*i;
    else if i>400
        Demand3(i)=200-0.5*(i-400);
        Demand10(i)=260-0.5*(i-400);
        Off_ramp(i)=330-0.9*(i-400);
    else 
        Demand3(i)=200;
        Demand10(i)=260;
        Off_ramp(i)=330;
        end
    end

    for k=1:14
        if k==3
            demand(k,i)=Demand3(i);
        else if k==8
                off_ramp(k,i)=Off_ramp(i);
            else if k==10
                    demand(k,i)=Demand10(i);
                else
                    demand(k,i)=0;
                    off_ramp(k,i)=0;
                end
            end
        end

         
        if ((density_d-density(3,i))>=0)
            On_ramp3(i)=fb_gain*(density_d-density(3,i));
        else On_ramp3(i)=ramp_min;
        end
        if (On_ramp3(i)<=(queue3(i)/T+Demand3(i)))                         %  if queue3(i)>=queue_max
                                                                           %     on_ramp(3,i)=ramp_max;
                                                                           %  else if (queue3(i)<queue_max)&&(on_ramp1(3,i)<=(queue3(i)/T+demand3(i)))
                on_ramp(3,i)=On_ramp3(i);                                  %        on_ramp(3,i)=on_ramp1(3,i);
                                                                           %       else
            else
                on_ramp(3,i)=queue3(i)/T+Demand3(i);                       %           on_ramp(3,i)=queue3(i)/T+demand3(i);
        end
        
        if ((density_d-density(10,i))>=0)
            On_ramp10(i)=fb_gain*(density_d-density(10,i));
        else On_ramp10(i)=ramp_min;
        end
        if (On_ramp10(i)<=(queue10(i)/T+Demand10(i)))
            on_ramp(10,i)=On_ramp10(i);
        else
            on_ramp(10,i)=queue10(i)/T+Demand10(k);
        end

        if k==1
            flow(k,i+1)=Flow;                                              %=Flow(i+1);here we need to define a mainstream flow Flow(i) that illustrate the peak hour flow
            velocity(k,i+1)=50;            
            density(k,i+1)=30;                                             %density(k,i+1)=flow(k,i+1)/velocity(k,i+i);             
        else if k==14
                density(k,i+1)=density(k-1,i);
                velocity(k,i+1)=velocity(k-1,i);
                flow(k,i+1)=flow(k-1,i);
            else 
                the_density(k,i)=velocity_free*(1-(density(k,i)/density_jam).^1.8).^1.7;
                
                flow(k,i)=alpha*density(k,i)*velocity(k,i)+(1-alpha)*(velocity(k+1,i)*density(k+1,i));
                
                density(k,i+1)=density(k,i)+(T/L)*(flow(k-1,i)-flow(k,i)+on_ramp(k,i)-off_ramp(k,i));
                
                velocity(k,i+1)=velocity(k,i)+(T/tao)*(the_density(k,i)-velocity(k,i))+(T/L)*velocity(k,i)*(velocity(k-1,i)-velocity(k,i))...
               -((gamma*T)/(L*tao))*(density(k+1,i)-density(k,i))/(density(k,i)+k_constant);
           
                queue3(i+1)=queue3(i)+T*(Demand3(i)-on_ramp(3,i));
                queue10(i+1)=queue10(i)+T*(Demand10(i)-on_ramp(10,i));
            end
        end
    end
end
figure;
i=2:1:13;mesh(density);xlabel('Time(15sec)');ylabel('Section');zlabel('Density');
figure;
i=1:1:600;
plot(density(3,i)*3,'k');hold on;
plot(queue3(i),'g');hold on;
plot(Demand3(i),'r');hold on;
plot(on_ramp(3,i),'c');hold on;
plot(density(10,i)*3,'k-.');hold on;
plot(queue10(i),'g-.');hold on;
plot(Demand10(i),'r-.');hold on;
plot(on_ramp(10,i),'c-.');hold on;
plot(off_ramp(8,i),'M');
xlabel('Time(15sec.)');