controller_mirm(densities, flows, celldata, ts, idx)

function new_flows = controller_mirm(densities, flows, celldata, ts, idx)
% CONTROLLER_MIRM - implementation of Modified "Ideal" Ramp Metering strategy,
%                   which takes into account downstream cell capacity and
%                   off-ramp flow in current cell.
%
% Call:   new_flows = controller_mirm(densities, flows, celldata, ts, idx)
%
% Parameters:
%             densities - vector of densities;
%             flows     - vector of on-ramp flows, same size as 'densities';
%             celldata  - array of freeway cell structures, whose length
%                         must be the same as size of 'densities';
%             ts        - sampling period;
%             idx       - index of the cell with our on-ramp.
%
% Returns:   new_flows - updated vector of on-ramp flows.
%
% Last modified:   10/20/2006.

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

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

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

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

if (idx < 1) | (idx > N)
  error('CONTROLLER_MIRM: invalid cell index.');
end

% initialize flow vector with current values
new_flows = flows;

% length of i-th cell
l         = abs(celldata(idx).PMend - celldata(idx).PMstart);
% given percentage of critical number of vehicles in i-th cell
nc        = celldata(idx).ORmlcontroller.ratio * celldata(idx).FDrhocrit * l;
% current number of vehicles in i-th cell
n         = densities(idx, 1) * l;
% on-ramp flow that i-th cell allows
f1        = (nc - n) / ts;

if (idx < N) & (celldata(idx+1).FDrhocrit >= densities(idx+1, 1))
  if celldata(idx).FRbeta < 1
    if celldata(idx).FRbeta < 1e-7
      F = celldata(idx+1).FDfmax;
    else
      b = (1 - celldata(idx).FRbeta) / celldata(idx).FRbeta;
      F = min([celldata(idx+1).FDfmax (b*celldata(idx).FRfmax)]);
    end
    b = 1 / (1 - celldata(idx).FRbeta);
    f2 = b * F  -  celldata(idx).FDfmax;
  else
    f2 = celldata(idx).FRfmax - celldata(idx).FDfmax;
  end
else
    f2 = 0;
end

f                 = max([0 f1 f2 (f1+f2)]);
f                 = min([f celldata(idx).ORmlcontroller.Cmax]);
f                 = max([f celldata(idx).ORmlcontroller.Cmin]);
new_flows(idx, 1) = f;  % update flow vector

return;