function rsys = remove_positive_numerator_time_constants(sys, tauc)
% REMOVE_POSITIVE_NUMERATOR_TIME_CONSTANTS  Cancels positive
% numerator time constants in LTI models against neighboring denominator
% terms in the model by applying the SIMC model reduction rules for such
% systems.
%
% RSYS = REMOVE_POSITIVE_NUMERATOR_TIME_CONSTANTS(SYS, TAUC) returns the
% simplified model, RSYS, where SYS is the original process model and TAUC
% is the tuning parameter (desired closed-loop time constant).
%
% See also REDUCE_PLANT_MODEL.

% Author(s): Bora Eryilmaz, The MathWorks, Inc.

[Z,P] = zpkdata(sys, 'v');

% Positive numerator time constants in decreasing order.
Zrealneg = Z(real(Z)<0 & ~imag(Z));
Tj0 = sort(-1./Zrealneg, 'descend');

% Lag time constants in descreasing order.
Prealneg = P(real(P)<0 & ~imag(P));
taui0 = sort(-1./Prealneg, 'descend');

% Remove positive numerator time constants.
rsys = zpk(sys);
for ct = 1:length(Tj0)
  T0 = Tj0(ct);

  % Find tau0a & tau0b such that tau0a >= T0 >= tau0b.
  tau0a = taui0(find(T0 <= taui0,1,'last'));
  tau0b = taui0(find(T0 >= taui0,1));

  % Find closest denominator time constant
  if isempty(tau0a) && isempty(tau0b)
    % No denominator term to cancel against.
    warning('SIMC:incomplete', ...
	    'Not all positive numerator time constants were removed.')
    break;
  elseif isempty(tau0a) || (~isempty(tau0b) && ...
			    ((T0 < sqrt(tau0a*tau0b)) && (T0 < 1.6*tau0b)))
    % T0/tau0b < tau0a/T0 and T0/tau0b < 1.6
    tau0 = tau0b;
  else
    tau0 = tau0a;
  end

  % Add approximation term and remove lead-lag term from the model.
  taut0 = [];
  if (T0 >= tau0) && (tau0 >= tauc)
    % Rule T1
    Gpz = zpk(-1/tau0,-1/T0,1);
  elseif (T0 >= tauc) && (tauc >= tau0)
    % Rule T1a
    Gpz = zpk(-1/tau0,-1/T0,tau0/tauc);
  elseif (tauc >= T0) && (T0 >= tau0)
    % Rule T1b
    Gpz = zpk(-1/tau0,-1/T0,tau0/T0);
  elseif (tau0 >= T0) && (T0 >= 5*tauc)
    % Rule T2
    Gpz = zpk(-1/tau0,-1/T0,1);
  elseif min(tau0,5*tauc) >= T0
    % Rule T3
    taut = min(tau0,5*tauc);
    Gpz = zpk(-1/tau0,[-1/T0 -1/(taut-T0)],taut/T0/(taut-T0));
    taut0 = taut-T0;
  else
    error('Unsupported case for a positive numerator time constant.')
  end

  % New model with positive time constant removed.
  rsys = minreal(rsys*Gpz);

  % Remove tau0 from the list of lag time constants, taui0, and
  % add the new time constant (from Rule T3) if needed.
  idx = find(tau0==taui0,1);
  taui0(idx) = [];
  taui0 = sort([taui0; taut0], 'descend');
end
end