function [K, CL, gamma, info] = h2syn_are(P, nmeas, ncont)
%h2syn_are Synthesizes the H2-norm minimizing dynamic output feedback controller for a regular generalized plant P.
%   P must have no (Cz, A)-unobservable or (A, Bw)-uncontrollable poles on the imaginary axis. Also P needs to be
%   (A, Bu)-stabilizable and  (Cy, A)-detectable.

arguments
    P lti
    nmeas (1,1) {mustBeInteger}
    ncont (1,1) {mustBeInteger}
end

%- Helper variables
nx = order(P);
[ny, nu] = size(P);  % Beware of conflict in defintion of nu,ny and later use.
widx = 1:(nu-ncont); nw = widx(end);
uidx = (widx(end)+1):nu; nu = uidx(end);
zidx = 1:(ny-nmeas); nz = zidx(end);
yidx = (zidx(end)+1):ny; ny = yidx(end);

A = P.A;
Bw = P.B(:,widx);  % Also known as B1
Bu = P.B(:,uidx);  % Also known as B
Cz = P.C(zidx,:);  % Also known as C1
Cy = P.C(yidx,:);  % Also known as C
Dzw = P.D(zidx, widx); assert(~any(Dzw, 'all'), "Non-zero feedthrough from z to w. Closed-loop has infinite H2-norm.")
Dzu = P.D(zidx, uidx); assert(rank(Dzu) >= min(size(Dzu)), "Dzu is not full column rank. Problem is not regular.");  % Also known as E or D1 or D12
Dyw = P.D(yidx, widx); assert(rank(Dyw) >= min(size(Dyw)), "Dyw is not full row rank. Problem is not regular.");  % Also known as F or D2 or D21
Dyu = P.D(yidx, uidx); assert(~any(Dyu, 'all'), "Non-zero feedthrough from y to u cannot be non-zero. Required transforms not implemented.")

Pyw = P(yidx, widx);
Pzu = P(zidx, uidx);

%- Solve
info = struct();
info.Pmax = icare(A, Bu, Cz' * Cz, Dzu'*Dzu, Cz'*Dzu, [], []);
assert(~isempty(info.Pmax), "ARE solving failed!");
info.Qmax = icare(A', Cy', Bw * Bw', Dyw*Dyw', Bw*Dyw', [], []);
assert(~isempty(info.Qmax), "Dual ARE solving failed!")
info.F = -(Dzu' * Dzu) \ (Dzu' * Cz + Bu' * info.Pmax);
info.E = -(info.Qmax * Cy' + Bw * Dyw') / (Dyw * Dyw');
%info.Jopt = sqrt(trace(Bw' * info.Pmax * Bw) + trace(Bu' * info.Pmax' * info.Qmax * info.Pmax * Bu));  % Only for standard form
info.Jopt = sqrt(trace(Bw' * info.Pmax * Bw) + trace((A' * info.Pmax' + info.Pmax' * A + Cz' * Cz) * info.Qmax));

K = ss(A + Bu * info.F + info.E * Cy, -info.E, info.F, 0);
K.TimeUnit = P.TimeUnit;

CL = lft(P, K);
gamma = norm(CL);
end