clear all;

%% Time
Time  = 660; %Time in minutes     
t = 0:Time;

%% Research on in Heuvelgalerie and internet
% People going in per minute
People_in_quiet                 = 1;
People_in_medium                = 1;
People_in_busy                  = 1;
People_in                       = People_in_quiet;    

Infected_percent                = 0.01;          % Percent of infected people
Time_measurement_station        = 1;            % Time people stay at the measurement station
Wait_time                       = 0;
%% Get optimum from matlab
Number_of_measurement_station   = 1;            % Max(5)
Random_percent                  = 0;            % Random people percentage that go to the measurement station per minute
Chance_infected_is_detected     = 0.8;          % Chance that a infected person is indeed picked out
Chance_healty_is_detected       = 0.02;         % Chance that a healthy person is detected

%% Define constants
m(1:5)                                  = 0;    % Measurement station is being used
w(1:5)                                  = 0;    % Time the measurement station is already used
People_measurement_station(1:Time+1)    = 0;    % total amount of people in the row for the measurement stations and inside of them
Row_length(1:Time+1)                    = 0;    % people waiting to use a measurement station
People(1:People_in)                     = 0;    % array of all the people, infected and healthy
Row = [];
D (1:Time+1)=0;
People_for_random                       = 2;
People_random = 0;
Person_random = 0;
%% Model
for Day_of_week = 1:7
    for n = 1:Time+1
        if (n==1)
            %% initial values at t = 0;
            People_measurement_station(n) = 0; 
            Row_length(n) = 0;
            m(1:5)                                  = 0;    % Measurement station is being used
            w(1:5)                                  = 0;    % Time the measurement station is already used
            Row                                     = [];
            People_in = People_in_quiet;
            Wait_time                               = 0;
        end
        
        %monday
        if (Day_of_week==1)&& (n == 181)
            People_in = People_in_medium;
        end
        if (Day_of_week==1)&& (n == 481)
            People_in = People_in_quiet;
        end
        if (Day_of_week==1)&& (n == 661)
            Mo_d = D;
            Mo_r = Row_length;
            Mo_p = People_measurement_station;
            Mo_Wait = Wait_time;
        end
        %Tuesday
        if (Day_of_week==2)&& (n == 121)
            People_in = People_in_medium;
        end
        if (Day_of_week==2)&& (n == 481)
            People_in = People_in_quiet;
        end
        if (Day_of_week==2)&& (n == 661)
            Tu_d = D;
            Tu_r = Row_length;
            Tu_p = People_measurement_station;
            Tu_Wait = Wait_time;
        end
        
        %Wednesday
        if (Day_of_week==3)&& (n == 241)
            People_in = People_in_medium;
        end
        if (Day_of_week==3)&& (n == 481)
            People_in = People_in_quiet;
        end
        if (Day_of_week==3)&& (n == 661)
            We_d = D;
            We_r = Row_length;
            We_p = People_measurement_station;
            We_Wait = Wait_time;
        end
        
        %Thursday
        if (Day_of_week==4)&& (n == 121)
            People_in = People_in_medium;
        end
        if (Day_of_week==4)&& (n == 481)
            People_in = People_in_quiet;
        end
        if (Day_of_week==4)&& (n == 661)
            Th_d = D;
            Th_r = Row_length;
            Th_p = People_measurement_station;
            Th_Wait = Wait_time;
        end
        %Friday
        if (Day_of_week==5)&& (n == 181)
            People_in = People_in_medium;
        end
        if (Day_of_week==5)&& (n == 661)
            Fr_d = D;
            Fr_r = Row_length;
            Fr_p = People_measurement_station;
            Fr_Wait = Wait_time;
        end
        
        %Saturday
        if (Day_of_week==6)&& (n == 121)
            People_in = People_in_medium;
        end
        if (Day_of_week==6)&& (n == 241)
            People_in = People_in_busy;
        end
        if (Day_of_week==6)&& (n == 481)
            People_in = People_in_medium;
        end
        if (Day_of_week==6)&& (n == 541)
            People_in = People_in_quiet;
        end
        if (Day_of_week==6)&& (n == 661)
            Sa_d = D;
            Sa_r = Row_length;
            Sa_p = People_measurement_station;
            Sa_Wait = Wait_time;
        end
        
        %Sunday
        if (Day_of_week==7)&& (n == 181)
            People_in = People_in_medium;
        end
        if (Day_of_week==7)&& (n == 241)
            People_in = People_in_busy;
        end
        if (Day_of_week==7)&& (n == 481)
            People_in = People_in_quiet;
        end
        if (Day_of_week==7)&& (n == 661)
            Su_d = D;
            Su_r = Row_length;
            Su_p = People_measurement_station;
            Su_Wait = Wait_time;
        end
        
        

        if (n>=2)
            %% Calculate infected and healthy people
            People(1:People_in) = 0;
            People = People(1:People_in);
            x = rand(1,People_in);
            for k = 1:People_in   %% people healthy or not
                if (x(k)>= Infected_percent)
                    People(k) = 0;                          % healthy people     
                else
                    People(k) = 1;                          % Infected people
                end     
            end

            People_infected = sum(People(:)==1); % Infected people per minute
            People_healthy = sum(People(:)==0);

            Detect_I(1:People_infected)                   = 0;
            Detect_H(1:People_healthy)                   = 0;
            Detect_I = Detect_I(1:People_infected);
            Detect_H = Detect_H(1:People_healthy);


    %% Calculate number of detected people

            x = rand(1,People_infected);
            for k = 1:People_infected
                if (x(k)>= Chance_infected_is_detected)
                    Detect_I(k) = 0;                         % infected people not detected     
                else
                    Detect_I(k) = 1;                         % Infected people detected
                end     
            end

            x = rand(1,People_healthy);
            for k = 1:People_healthy
                if (x(k)>= Chance_healty_is_detected)
                Detect_H(k) = 0;                         % healthy people not detected     
                else
                    Detect_H(k) = 1;                         % healthy people detected
                end     
            end

            People_detected = sum(Detect_I(:)==1)+sum(Detect_H(:)==1); % Detected people per minute

            D(n-1) = People_detected;
   
            

    %% declare people being tested in the measurement station and the row
            if (Person_random >= 2 && People_detected == 0)
                People_random = 1;
                Person_random = Person_random -2;
            end
            People_measurement_station(n) = People_detected + People_measurement_station(n-1); 
            Row_length(n)= People_detected + Row_length(n-1);
            
            People_random = 0;
            
            
            
            x = randperm(People_detected, sum(Detect_I(:)==1));
            Row_add(1:People_detected) = 0;
            Row_add = Row_add(1:People_detected);
            Row_add(x(:)) = 1;
            Row = [Row, Row_add];
            
            %% people finished the checkup
            % Measurement station 1
            if  (w(1) == Time_measurement_station) && (Number_of_measurement_station >= 1)
                m(1) = 0;
                w(1) = 0;
                People_measurement_station(n) = People_measurement_station(n)-1;
            end
            % Measurement station 2
            if  (w(2) == Time_measurement_station) && (Number_of_measurement_station >= 2)
                m(2) = 0;
                w(2) = 0;
                People_measurement_station(n) = People_measurement_station(n)-1;
            end
            % Measurement station 3
            if  (w(3) == Time_measurement_station) && (Number_of_measurement_station >= 3)
                m(3) = 0;
                w(3) = 0;
                People_measurement_station(n) = People_measurement_station(n)-1;
            end
            % Measurement station 4
            if  (w(4) == Time_measurement_station) && (Number_of_measurement_station >= 4)
                m(4) = 0;
                w(4) = 0;
                People_measurement_station(n) = People_measurement_station(n)-1;
            end
            % Measurement station 5
            if  (w(5) == Time_measurement_station) && (Number_of_measurement_station >= 5)
                m(5) = 0;
                w(5) = 0;
                People_measurement_station(n) = People_measurement_station(n)-1;
            end
            %% if there is a person in the row and a medical station is free person goes inside medical station
            % Measurement station 1
            if  ((Row_length(n) >= 1) && (m(1) < 1)) &&(Number_of_measurement_station >= 1)
                m(1) = 1;
                Row_length(n) = Row_length(n)-1;
                Row = Row(2:end);
            end
            % Measurement station 2            
            if  ((Row_length(n) >= 1) && (m(2) < 1))&& (Number_of_measurement_station >= 2)
                m(2) = 1;
                Row_length(n) = Row_length(n)-1;
                Row = Row(2:end);
            end 
            % Measurement station 3  
            if  ((Row_length(n) >= 1) && (m(3) < 1))&& (Number_of_measurement_station >= 3)
                m(3) = 1;
                Row_length(n) = Row_length(n)-1;
                Row = Row(2:end);
            end
            % Measurement station 4
            if  ((Row_length(n) >= 1) && (m(4) < 1))&& (Number_of_measurement_station >= 4)
                m(4) = 1;
                Row_length(n) = Row_length(n)-1;
                Row = Row(2:end);
            end
            % Measurement station 5
            if  ((Row_length(n) >= 1) && (m(5) < 1))&& (Number_of_measurement_station >= 5)
                m(5) = 1;
                Row_length(n) = Row_length(n)-1;
                Row = Row(2:end);
            end

            %% Waiting time increase
            % Measurement station 1
            if (m(1)==1)
                w(1)=w(1)+1;
            end
            % Measurement station 2
            if (m(2)==1)
                w(2)=w(2)+1;
            end
            % Measurement station 3
            if (m(3)==1)
                w(3)=w(3)+1;
            end
            % Measurement station 4
            if (m(4)==1)
                w(4)=w(4)+1;
            end
            % Measurement station 5
            if (m(5)==1)
                w(5)=w(5)+1;
            end
            if (Row_length(n) == 0)
                Person_random = Person_random + People_detected;
            end
            Wait_time(n) = (Row_length(n) * Time_measurement_station - w(1)-w(2)-w(3)-w(4)-w(5))/Number_of_measurement_station;
        end
    end
end


% Monday
subplot(7,1,1)
stairs(t/60+10,Mo_d)
hold on
plot(t/60+10,Mo_r)
plot(t/60+10,Mo_p)
hold off
legend('detected','Row','People total');

% Tuesday
subplot(7,1,2)
stairs(t/60+10,Tu_d)
hold on
plot(t/60+10,Tu_r)
plot(t/60+10,Tu_p)
hold off
legend('detected','Row','People total');

%Wednesday
subplot(7,1,3)
stairs(t/60+10,We_d)
hold on
plot(t/60+10,We_r)
plot(t/60+10,We_p)
hold off
legend('detected','Row','People total');

%Thursday
subplot(7,1,4)
stairs(t/60+10,Th_d)
hold on
plot(t/60+10,Th_r)
plot(t/60+10,Th_p)
hold off
legend('detected','Row','People total');

%Friday
subplot(7,1,5)
stairs(t/60+10,Fr_d)
hold on
plot(t/60+10,Fr_r)
plot(t/60+10,Fr_p)
hold off
legend('detected','Row','People total');

%Saturday
subplot(7,1,6)
stairs(t/60+10,Sa_d)
hold on
plot(t/60+10,Sa_r)
plot(t/60+10,Sa_p)
hold off
legend('detected','Row','People total');

%Saturday
subplot(7,1,6)
stairs(t,Sa_d)
hold on
plot(t,Sa_r)
plot(t,Sa_p)
hold off
legend('detected','Row','People total');

%Sunday
subplot(7,1,7)
stairs(t/60+10,Su_d)
hold on
plot(t/60+10,Su_r)
plot(t/60+10,Su_p)
hold off
legend('detected','Row','People total');

