function chk=h2synchk(P,K)
% function chk=h2synchk(P,K)
%
% checking K for H2-optimality with plant P
[A,B,C,D]=ssdata(P);
[Af,Bf,Cf,Df]=ssdata(K);
[mdf,ndf]=size(Df);
na=size(A,1);
naf=size(Af,1);
[md,nd]=size(D);
ndist=nd-mdf;
ncost=md-ndf;
if ndist<=0, error('K has too many inputs'); end
if ncost<=0, error('K has too many outputs'); end
X=[Af Bf;Cf Df];
a1=[zeros(na,naf) B(:,ndist+1:nd); eye(naf) zeros(naf,mdf)];
a2=[zeros(naf,na) eye(naf); C(ncost+1:md,:) zeros(ndf,naf)];
a0=[A zeros(na,naf);zeros(naf,na+naf)]+a1*X*a2;
mrg=max(abs(eig(a0)));      % stability margin
if mrg>0.999999, chk=Inf; return; end

b1=a1;
b2=[zeros(naf,ndist);D(ncost+1:nd,1:ndist)];
b0=[B(:,1:ndist);zeros(naf,ndist)]+b1*X*b2;
c1=[zeros(ncost,naf) D(1:ncost,ndist+1:nd)];
c2=a2;
c0=[C(1:ncost,:) zeros(ncost,naf)]+c1*X*c2;
d1=c1;
d2=b2;
d0=D(1:ncost,1:ndist)+d1*X*d2;

P=dlyap(a0,b0*b0');
Q=dlyap(a0',c0'*c0);
chk=norm(d2*d0'*d1+c2*P*c0'*c1+b2*b0'*Q*b1+a2*P*a0'*Q*a1);