compute_flows(densities, celldata, inflow, ts)

function flows = compute_flows(densities, celldata, inflow, ts)
% COMPUTE_FLOWS - computes and returns flows based on densities,
%                 split ratios and fundamental diagrams.
%
% Call:   flows = compute_flows(densities, celldata, inflow, ts)
%
% Parameters:
%             densities - vector of 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/16/2006.

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

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

if M ~= N
  error('COMPUTE_FLOWS: number of densities 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;
    % off-ramp split ratio in (i-1)th cell
    b    = celldata(i-1).FRbeta;
    % max flow that can be generated by (i-1)th cell
    val3 = (1 - b) * (v/l) * (n + g*r);
  end

  fl = min([val1 val2 val3]);

  % consistency check: main line flow cannot exceed ((1-b)/b)*FRfmax
  if i > 1
    b = celldata(i-1).FRbeta;
    if (b > 0) & (b < 1)
      if ((b/(1-b)) * fl) > celldata(i-1).FRfmax
        fl = ((1-b)/b) * celldata(i-1).FRfmax; 
      end
    elseif b > 0
      fl = 0;
    end
  end

  flows = [flows; fl];
end

return;