function [K,CL,GAM]=h2syn6245(p,nmeas,ncon)
% function [K,CL,GAM]=h2syn6245(p,nmeas,ncon)
% 
% 6.245 version of h2syn.m for DT: 
% H2 optimal controller K=ss(Af,Bf,Cf,Df,-1)
% for open loop plant   p=ss(A,[B1 B2],[C1;C2],[D11 D12;D21 zeros(k,m)],-1)

[A,B,C,D]=ssdata(p);
nstate=size(A,1);
[nout,nin]=size(D);
ncost=nout-nmeas; 
if ncost<=0, error('number of measurements is too large'); end
ndist=nin-ncon;
if ndist<=0, error('number of controls is too large'); end
C1=C(1:ncost,:);
C2=C(ncost+1:nout,:);
B1=B(:,1:ndist);
B2=B(:,ndist+1:nin);
D11=D(1:ncost,1:ndist);
D12=D(1:ncost,ndist+1:nin);
D21=D(ncost+1:nout,1:ndist);

[Pc,~,F,rep] = dare(A,B2,C1'*C1,D12'*D12,C1'*D12);
if rep==-1, error('ker[A-zI B2;C1 D12] not {0} for some |z|=1'); end
if rep==-2, error('range[A-zI B2] not full for some |z|>=1'); end
F=-F;
[Po,~,L,rep] = dare(A',C2',B1*B1',D21*D21',B1*D21');
if rep==-1, error('range[A-zI B1;C2 D21] not full for some |z|=1'); end
if rep==-2, error('ker[A-zI;C2] not {0} for some |z|>=1'); end
L=-L';
Ec=chol6245(Pc); 
Eo=chol6245(Po);
Gc=chol6245(D12'*D12+B2'*Pc*B2);

% reshape(a*x,size(a,1)*size(x,2),1)=kron(eye(size(x,2)),a)*x(:)
%                                   =kron(x.',eye(size(a,1)))*a(:)
% form matrix of the map 
%  K -> (Ec*B2*K*D21,D12*K*D21,Gc*K*C2*Eo'), 
% comp. to (-Ec*B1,-D11,Gc*F*Eo')
M1=kron(eye(ndist),Ec*B2)*kron(D21',eye(ncon)); 
y1=-reshape(Ec*B1,nstate*ndist,1);
M2=kron(eye(ndist),D12)*kron(D21',eye(ncon));
M3=kron(eye(nstate),Gc)*kron(Eo*C2',eye(ncon));
y3=reshape(Gc*F*Eo',nstate*ncon,1);
Df = reshape([M1;M2;M3]\[y1;-D11(:);y3],ncon,nmeas);
Af=A+B2*F+L*C2-B2*Df*C2;
Bf=B2*Df-L;
Cf=F-Df*C2;
K=ss(Af,Bf,Cf,Df,-1);
Acl=[A+B2*Df*C2 B2*Cf;Bf*C2 Af];
Bcl=[B1+B2*Df*D21;Bf*D21];
Ccl=[C1+D12*Df*C2, D12*Cf];
Dcl=D11+D12*Df*D21;
CL=ss(Acl,Bcl,Ccl,Dcl,-1);
GAM=norm([norm(Ec*(B1+B2*Df*D21),'fro'); ...
    norm(Gc*(F-Df*C2)*Eo','fro'); ...
    norm(D11+D12*Df*D21,'fro')]);