compute_densities(densities, flows, frflows, celldata, outflow, ts)

function newden = compute_densities(densities, flows, frflows, celldata, outflow, ts)
% COMPUTE_DENSITIES - implementation of conservation equation.
%
% Call:   newden = compute_densities(densities, flows, frflows, celldata, outflow, ts)
%
% Parameters:
%             densities - vector of densities;
%             flows     - vector of flows, must have the same size as 'densities';
%             frflows   - vector of off-ramp flows, must have the same size as 'densities';
%             celldata  - array of cell data structures,
%                         whose length must be the same as size of 'densities';
%             outflow   - how much can be allowed through the right boundary;
%             ts        - sampling period.
%
% Returns:   newden - vector of densities for time step (k + 1).
%
% Last modified:   09/14/2006.

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

K = size(densities, 1);
L = size(flows, 1);
M = size(frflows, 1);
N = size(celldata, 2);  % row vector

if K ~= L
  error('COMPUTE_DENSITIES: number of densities does not match number of flows.');
end

if K ~= M
  error('COMPUTE_DENSITIES: number of densities does not match number of off-ramp flows.');
end

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

newden = [];

for i = 1:N
  l   = abs(celldata(i).PMend - celldata(i).PMstart);  % cell length
  d   = densities(i, 1);  % current density for i-th cell
  qin = flows(i, 1)  +  celldata(i).ORflow;  % total flow entering i-th cell

  if i == N  % last cell
    qout = get_rbflow(densities(N, 1), flows(N, 1), celldata(N), outflow, ts);
  else
    qout = flows(i+1, 1) + frflows(i, 1);
  end

  nd = d + (ts/l) * (qin - qout);  % new density for the i-th cell

  if nd >= celldata(i).FDrhojam
    nd = celldata(i).FDrhojam;
  elseif nd <= 0
    nd = 0;
  end

  newden = [newden; nd];

end


return;