%[Nl, Ml] = leftcopr(G, tol);
%
%  Finds the left normalized coprime fraction (Nl,Ml) of G = Ml^{-1}*Nl.
%  See book of Skogestad and Postlethwaite Chapter 4, page 118.
%
%  Written by: Kjetil  Havre, 7/3-1996.
%
function [Nl, Ml, Z, R] = leftcopr(G, tol)

   Nl=[]; Ml=[]; Z=[]; R=[];
   if nargin < 2
      tol = 1e-15;
   end 

   [a,b,c,d] = unpck(G);
   
   R = eye(size(d*d.')) + d*d.';
   S = eye(size(d.'*d)) + d.'*d;
   
   A1 = (a - b*pinv(S)*d.'*c); 
   M  = -c.'*pinv(R)*c;
   Q  = -b*pinv(S)*b.';
   
   H1 = [A1' M; Q -A1];
   
   %MREH1 = min(abs(real(eig(H1))))
  
   [Z1, Z2, FAIL, reig_min] = ric_schr(H1, tol);
   %FAIL = FAIL
   %reig_min = reig_min
   if FAIL == 1
      disp('Eigenvalues on jw-axis.');
      return;
   end
   
   if FAIL == 2
      disp('Unequal number of stable and unstable eigenvalues.');
      return;  
   end
   
   if FAIL == 2
      disp('Both unequal number of stable and unstable eigenvalues and');
      disp('eigenvalues on jw-axis.');
      return;     
   end
   
   Z = Z2/Z1;
   H = -(b*d.' + Z*c.')*pinv(R);
   R2 = pinv(sqrt(R));
   
   A = a+H*c; Bn = b+H*d;; Bm = H; C = R2*c; Dn = R2*d; Dm = R2;   
   Nl = pck(A, Bn, C, Dn); Ml = pck(A, Bm, C, Dm);
   
end
