function [N,M]=ncopfac(a,b,c,d)
%NCOPFAC    [M,N]=ncopfac(A,B,C,D) Uses the Mu toolbox to find
%           normalized coprime factors of system [A,B,C,D] using (4.26).
% 11Nov99: Corrected by replacing sqrt by sqrtm
%           
S=eye(size(d'*d))+d'*d; 
R=eye(size(d*d'))+d*d';
A1 = a-b*inv(S)*d'*c;
R1 = c'*inv(R)*c;
Q1 = b*inv(S)*b';
[z1,z2,fail,reig_min] = ric_schr([A1' -R1; -Q1 -A1]);
Z = z2/z1;

% Alternative: aresolv in Robust control toolbox:
% [z1,z2,eig,zerr,zwellposed,Z] = aresolv(A1',Q1,R1);

H = -(b*d' + Z*c')*inv(R);
A = a + H*c;
Bn = b + H*d;
Bm = H;
C = inv(sqrtm(R))*c;
Dn = inv(sqrtm(R))*d;
Dm = inv(sqrtm(R));
N = pck(A,Bn,C,Dn);
M = pck(A,Bm,C,Dm);
