% Deterministic queueing simulation of an airport facility
% Treats demand and capacity as continuous functions
% (time dependent)
% 
% Integrates the demand and supply flow rates to get
% cumulative supply and demand curves
%
% Basic state equations =>

% p(1) = L(t) = L(t-1) + [demand(t) - supply(t)] (dt)           % queue length 
% p(2) = A(t) = A(t-1) + [L(t)] (dt)                                      % area under L(t) curve
%
%  Demand / Supply
%		| 
%		|________			Demand 	
%		|   	 \			__________
%		|_________\_______/___________\______ Supply	
%       |			-----			   \
%		|		   
%		|		    
%		|___________________________________
%			    			Time
%
%
% Programmer: Toni Trani
% Date: June 2013
% 
% Calls: fqueue_deterministic_mod.m
%
%  fqueue_deterministic_mod.m computes the integrals of the demand and supply rates
%  of change

clear
%
% Define global variables

global demand capacity time

% Enter demand function as an array of values over time

% time = time vector containing discrete points in time
% demand = vector with values of demand rate  (entities/time) observed in the system
% capacity = vector with values of capacity rate (entities/time) for the system

	time = [0 1 1.00001 2 2.00001 3 4 ];                           % values of time (time vector)
	demand = [1500 1500 500 500 500 500 500];             % values of demand over time
	capacity = [1000 1000 1000 1000 1500 1500 1500 ];  % values of capacity over time

% Compute min and maximum values for proper scaling in plots
mintime 	= min(time);
maxtime 	= max(time+3);
maxd	 	= max(demand);
maxc 		= max(capacity);
mind	 	= min(demand);
minc 		= min(capacity);
scale		= round(.2 *(maxc+maxd)/2)
minplot = min(minc,mind) - scale;
maxplot	= max(maxc,maxd) + scale;

po = [0 0 0 0 0];			% intial number of vehicles in the queue (assumed to be zero)
to = mintime;
tf = maxtime;
tspan = [to tf];

% where:
% to is the initial time to solve this equation
% tf is the final time
% tspan is the time span to solve the simulation

% Call the ODE solver in Matlab and call function fqueue_deterministic that
% contains the rates of change of demand and capacity

[t,p] = ode23('fqueue_deterministic_mod',tspan,po);

% Compute statistics

Ltmax 		= max(p(:,1));
tdelay 		= max(p(:,2));

  clc
  disp([blanks(5),'Deterministic Queueing Model '])
  disp(' ')
  disp(' ')
  disp([blanks(5),'  Simulation Period (hours) = ', num2str(maxtime)])

  disp(' ')
  
  disp(' ')
  disp([blanks(5),'  Total delay (aircraft-hour) = ', num2str(tdelay)])
  disp([blanks(5),'  Max queue length (aircraft) = ', num2str(Ltmax)])
  disp(' ')


% Plot the demand and supply functions

figure

subplot(2,1,1)
plot(time,demand,'b',time,capacity,'--k')
xlabel('Time (hours)')
ylabel('Demand or Capacity Rates (Entities/hour)')
axis([mintime maxtime minplot maxplot])
legend('Demand Rate','Capacity Rate')
grid

subplot(2,1,2)
plot(t,p(:,1),'b')
xlabel('Time (Hours)')
ylabel('Entity in Queue')
grid

% pause

figure

% Plot the results of the numerical integration procedure

subplot(2,1,1)
plot(t,p(:,1),'b')
xlabel('Time (hours)')
ylabel('Entities in Queue')
grid

subplot(2,1,2)
plot(t,p(:,2),'k')
xlabel('Time (hours)')
ylabel('Total Delay (Entities-hour)')
grid


% new figure to plot integral or demand and supplu independently
figure

% Plot the results of the numerical integration procedure for demand
% function

plot(t,p(:,4),'o-k',t,p(:,5),'^--r')
xlabel('Time (hours)')
ylabel('Cumulative Number of Entities ')
grid


