compute_flows_2(densities, frflows, celldata, inflow, ts)

function flows = compute_flows_2(densities, frflows, celldata, inflow, ts)
% COMPUTE_FLOWS_2 - computes and returns flows based on densities,
%                   off-ramp flows and fundamental diagrams.
%
% Call:   flows = compute_flows_2(densities, frflows, celldata, inflow, ts)
%
% Parameters:
%             densities - vector of densities;
%             frflows   - vector of off-ramp flows, same size as 'densities';
%             celldata  - array of freeway cell structures, whose length
%                         must be the same as size of 'densities';
%             inflow    - demand on the left boundary of the freeway section;
%             ts        - sampling time.
%
% Returns:   flows - vector of flows.
%
% Last modified:   09/29/2006.

%
% Alex Kurzhanskiy   <[e-mail]>
%

L = size(densities, 1);
M = size(frflows, 1);
N = size(celldata, 2);

if L ~= N
  error('COMPUTE_FLOWS_2: number of densities does not match number of cells.');
end

if M ~= N
  error('COMPUTE_FLOWS_2: number of off-ramp flows does not match number of cells.');
end

flows = [];

for i = 1:N
  % i-th cell length
  l    = abs(celldata(i).PMend - celldata(i).PMstart);
  % number of vehicles that creates jam in i-th cell
  nJ   = celldata(i).FDrhojam * l;
  % actual number of vehicles in i-th cell
  n    = densities(i, 1) * l;
  % number of vehicles coming from on-ramp into i-th cell
  r    = celldata(i).ORflow * ts;
  % on-ramp flow blending coefficient in i-th cell
  g    = celldata(i).ORgamma;
  % congestion wave speed in i-th cell
  w    = celldata(i).FDfmax / (celldata(i).FDrhojam - celldata(i).FDrhocrit);
  % max flow that i-th cell allows to enter when congested
  val2 = max([0 ((w/l) * (nJ - (n + g*r)))]);

  if i == 1
    % max flow that can enter i-th cell 
    val1 = celldata(i).FDfmax;
    % demand on the left boundary of freeway section
    val3 = inflow;
  else
    % max flow that can enter i-th cell 
    val1 = celldata(i-1).FDfmax;
    % (i-1)th cell length
    l    = abs(celldata(i-1).PMend - celldata(i-1).PMstart);
    % actual number of vehicles in (i-1)th cell
    n    = densities(i-1, 1) * l;
    % number of vehicles coming from on-ramp into (i-1)th cell
    r    = celldata(i-1).ORflow * ts;
    % on-ramp flow blending coefficient in (i-1)th cell
    g    = celldata(i-1).ORgamma;
    % free-flow speed in (i-1)th cell
    v    = celldata(i-1).FDfmax / celldata(i-1).FDrhocrit;
    % max flow that can be generated by (i-1)th cell
    val3 = max([0 ((v/l) * (n + g*r)  -  frflows(i-1, 1))]);
  end

  flows = [flows; min([val1 val2 val3])];
end

return;