global Q_LIMIT mean_interarrival mean_service mean_delays_required
global next_event_type num_custs_delayed num_delays_required num_events num_in_q server_status
global area_num_in_q area_server_status mean_interarrival mean_service time time_arrival time_last_event time_next_event total_of_delays

% Input Parameters
Q_LIMIT = 200;
mean_interarrival = 1.0;
mean_service = 0.5;
mean_delays_required = 100;
area_server_status=0; num_events = 2
INIT; % calls initialization routine

% run the simulation while more delays are still required
while(num_custs_delayed < mean_delays_required)
    TIMING;
    UPDATE;
    if(next_event_type==1)
        ARRIVE; % branch to arrival routine (function)
    elseif(next_event_type==2)
        DEPART; % branch to departure routine (function)
    end
end

REPORT % call report generator function
function INIT
global Q_LIMIT mean_interarrival mean_service mean_delays_required
global next_event_type num_custs_delayed num_delays_required num_events num_in_q server_status
global area_num_in_q area_server_status mean_interarrival mean_service time time_arrival time_last_event time_next_event total_of_delays
time=0.0;
server_status=0;
num_in_q = 0;
time_last_event = 0.0;
num_custs_delayed = 0;
total_of_delays = 0.0;
area_num_in_q = 0.0;
area_erve_status = 0.0;
time_next_event(1) = time+EXPON(mean_interarrival)
time_next_event(2) = 1.0 * exp(30) %disp(['init'])
end

function ARRIVE
global Q_LIMIT mean_interarrival mean_service mean_delays_required
global next_event_type num_custs_delayed num_delays_required num_events num_in_q server_status
global area_num_in_q area_server_status mean_interarrival mean_service time time_arrival time_last_event time_next_event total_of_delays
time_next_event(1) = time + EXPON(mean_interarrival);
% time
% check to see whether server is busy
if(server_status == 1) % server is busy
num_in_q = num_in_q + 1; % increase the no. of customers in queue
if (num_in_q) > Q_LIMIT
    disp(['num_in_q = ', num2str(num_in_q)]);
    disp(['Overflow of the array of time_arrival at ', num2str(time)]);
    pause
end
time_arrival(num_in_q) = time;
else % server is idle
    delay = 0.0;
    total_of_delays = total_of_delays + delay;
    num_custs_delayed = num_custs_delayed + 1;
    server_status = 1;
    time_next_event(2) = time + EXPON(mean_service);
end
end

function DEPART
global Q_LIMIT mean_interarrival mean_service mean_delays_required
global next_event_type num_custs_delayed num_delays_required num_events num_in_q server_status
global area_num_in_q area_server_status mean_interaarival mean_service time time_arrival time_last_event time_next_event total_of_delays
if(num_in_q == 0) % queue is empty. Make the server idle and eliminate the departure event from consideration
    server_status = 0; % idle
    time_next_event(2) = 1.0 * exp(30);
else % queue is not empty, so decrease the no. of customers in queue
    num_in_q = num_in_q - 1;
    delay = time - time_arrival(1);
    total_of_delays = total_of_delays + delay;
    num_custs_delayed = num_custs_delayed + 1;  %???????
    time_next_event(2) = time + EXPON(mean_service);
    % move each customers's arrival time in queue up one place
    for i = 1:num_in_q
        time_arrival(i)=time_arrival(i+1);
    end
end
end

function e=EXPON(mean)
u = rand(1);
e = - mean * log(u);
% uses the inverse transformation method to generate negativeexponential random numbers
end

function TIMING
global Q_LIMIT mean_interarrival mean_service mean_delays_required
global next_event_type num_custs_delayed num_delays_required num_events num_in_q server_status
global area_num_in_q area_server_status mean_service time time_arrival time_last_event time_next_event total_of_delays
min_time_next_event = 1.0*exp(29);
next_event_type = 3; % determine the event type of the next event to occur
for i=1:num_events
    if(time_next_event(i) < min_time_next_event)
        min_time_next_event = time_next_event(i);
        next_event_type = i;
    end;
end
if(next_event_type == 3)
    disp(['Event List Empty at Time ', num2str(time)]);
end
time = min_time_next_event;
end

function UPDATE
% Update area accumulated for time-average statistics
global Q_LIMIT mean_interarrival mean_service mean_delays_required
global next_event_type num_custs_delayed num_delays_required num_events num_in_q server_status
global area_num_in_q area_server_status mean_interarrival mean_service time time_arrival time_last_event time_next_event total_of_delays
time_since_last_event = time - time_last_event;
time_last_event = time;
area_num_in_q = area_num_in_q + num_in_q * time_since_last_event;
area_server_status = area_server_status + server_status * time_since_last_event;
end

function REPORT
global Q_LIMIT mean_interarrival mean_service mean_delays_required
global next_event_type num_custs_delayed num_delays_required num_events num_in_q server_status
global area_num_in_q area_server_status mean_interarrival mean_service time time_arrival time_last_event time_next_event total_of_delays
disp(['Average Delay in Queue : ', num2str(total_of_delays/num_custs_delayed),' minutes']);
disp(['Average number in Queue : ', num2str(area_num_in_q/time)])
disp(['Service Utilization : ', num2str(area_server_status/time)]);
disp(['Time Simulation Ended : ', num2str(time)]);
end