function [H,H0,INFO] = modsepREG(G,N,RegionFcn,opt)
% Region-based modal decomposition.

%   Author(s): P. Gahinet
%   Copyright 1986-2023 The MathWorks, Inc.
[ny,nu] = iosize(G);
if hasInternalDelay(G)
   throw(ltipack.utNoDelaySupport('modalsep',G.Ts,'internal'))
end

% Extract data
a = G.a;  b = G.b;  c = G.c;  d = G.d;  e = G.e;   Ts = G.Ts;
dx0 = G.dx0; x0 = G.x0;  u0 = G.u0;  y0 = G.y0;
nx = size(a,1);
EXPLICIT = isempty(e);
REAL = isreal(a) && isreal(e);

% Scale A,B,C to minimize loss of accuracy in Schur/QZ decomposition
% and reordering
if G.Scaled
   % No scaling necessary
   sL = ones(nx,1);  sR = sL;
else
   [a,b,c,e,sR,sL] = ltipack.xscale(a,b,c,d,e,Ts);
   dx0 = sL.*dx0;  x0 = sR.\x0;
end

% Schur/QZ decomposition
if EXPLICIT
   [u,a] = schur(a);
   p = ordeig(a);
else
   if REAL
      [a,e,q,z] = qz(a,e,'real');
   else
      [a,e,q,z] = qz(a,e);
   end
   p = ordeig(a,e);
end

% Assign region index to each eigenvalue
try
   clusters = zeros(nx,1);
   for ct=1:nx
      clusters(ct) = round(RegionFcn(p(ct)));
   end
catch  %#ok<CTCH>
   error(message('Control:transformation:modsep1'))
end
if any(clusters<1 | clusters>N)
   error(message('Control:transformation:modsep2'))
end

% Go complex if (A,E) real but different regions are assigned to conjugate pairs
if REAL
   idxp = find(imag(p)>0);
   [~,ia,ib] = intersect(p,conj(p(idxp)));
   if any(clusters(ia)~=clusters(idxp(ib)))
      if EXPLICIT
         [u,a] = rsf2csf(u,a);
      else
         [a,e,qq,zz] = qz(a,e);
         q = qq*q;  z = z*zz;
      end
   end
end

% Sort eigenvalues by increasing region index
try
   if EXPLICIT
      [u,a] = ordschur(u,a,-clusters);
      b = u'*b;  c = c*u;  
      dx0 = u'*dx0;  x0 = u'*x0;
   else
      [a,e,q,z] = ordqz(a,e,q,z,-clusters);
      b = q*b;  c = c*z;
      dx0 = q*dx0;  x0 = z'*x0;
   end
catch %#ok<CTCH>
   error(message('Control:transformation:modsep4'))
end
clusters = sort(clusters);
if isempty(clusters)
   blksize = [];
else
   blksize = diff([0;find(diff(clusters));numel(clusters)]);
end

% Block diagonalize
try
   if EXPLICIT
      [tR,a] = bdschur(a,[],blksize);
      aux = tR\[b dx0 x0];
      b = aux(:,1:nu);  dx0 = aux(:,nu+1);  x0 = aux(:,nu+2);
      c = c*tR;
   else
      [a,e,tL,tR] = bdqz(a,e,[],blksize);
      b = tL'*b;  dx0 = tL'*dx0;
      c = c*tR;
      if any(x0)
         x0 = tR\x0;
      end
   end
catch ME
   % Should not happen
   throw(ME)
end

% Initialize outputs
H0 = ltipack.ssdata([],zeros(0,nu),zeros(ny,0),d,[],Ts);
H0.u0 = u0;  H0.y0 = y0;  
H0.Delay = G.Delay;
% N-by-1 array H
h = ltipack.ssdata([],zeros(0,nu),zeros(ny,0),zeros(ny,nu),[],Ts);
h.u0 = u0;  h.y0 = y0;  
h.Delay = G.Delay;
H = repmat(h,[N 1]);

% Store factors
is = 0;
for j=1:numel(blksize)
   % Store j-th component
   bs = blksize(j);
   ix = is+1:is+bs;
   jsys = clusters(is+1);
   H(jsys).a = a(ix,ix);
   H(jsys).b = b(ix,:);
   H(jsys).c = c(:,ix);
   if ~EXPLICIT
      H(jsys).e = e(ix,ix);
   end
   H(j).dx0 = dx0(ix,:);
   H(j).x0 = x0(ix,:);
   is = is+bs;
end

% Construct block diagonalization state transformation T and its inverse Ti
% RE: The system with decoupled dynamics is (TL*A*TR,TL*B,C*TR,D)
if nargout>2
   INFO = struct('Mode',[],'DCGain',zeros(N,1),'TL',[],'TR',[]);
   if EXPLICIT
      INFO.TL = tR \ (sL .* u)';
      INFO.TR = (sR .* u) * tR;
   else
      INFO.TL = tL' * (q .* sL');
      INFO.TR = (sR .* z) * tR;
   end
   if Ts==0
      s = opt.DCFrequency;
   else
      s = exp(opt.DCFrequency*abs(Ts));
   end
   DCG = norm(evalfr(G,s));
   for j=1:N
      INFO.DCGain(j) = norm(evalfr(H(j),s))/DCG;
   end
end
