function [M0,W1,W2,Z] = getSectorData(Q,dQ)
% Computes factorizations for sector index computations
%    dQ=[] -> M0 = [],  Q = Z(s)' (W1*W1'-W2*W2') * Z(s)
%    Q,dQ  -> M0 = Q,  dQ = W2*W2',  W1=[],  Z=[]

%   Copyright 1986-2015 The MathWorks, Inc.
if isnumeric(Q) && ~(allfinite(Q) && isequal(Q,Q'))
   error(message('Control:analysis:getSectorIndex3'))
end
Z = [];
if isempty(dQ)
   % R-index
   if isa(Q,'numlti')
      % Factor Q(s) = Z(s)' * Q * Z(s) with Z(s) bi-stable
      try
         Q = ss(Q);
      catch
         error(message('Control:analysis:getSectorIndex12'))
      end
      try
         [Z,Q] = spectralfact(Q);
      catch ME
         error(message('Control:analysis:getSectorIndex14',ME.message))
      end
   end
   % Factorize Q = W1*W1' - W2*W2'
   [U,D] = schur(Q);
   D = diag(D);
   i1 = find(D>=0);
   i2 = find(D<0);
   M0 = [];
   W1 = lrscale(U(:,i1),[],sqrt(D(i1)));
   W2 = lrscale(U(:,i2),[],sqrt(-D(i2)));
else
   % Directional index for direction dQ
   if ~isnumeric(Q)
      % Directional index not supported for frequency-dependent Q
      error(message('Control:analysis:getSectorIndex5'))
   elseif ~(isnumeric(dQ) && isequal(size(dQ),size(Q)) && isequal(dQ,dQ'))
      error(message('Control:analysis:getSectorIndex4',size(Q,1)))
   end
   % Factorize dQ = W*W'
   [U,D] = schur(dQ);
   D = diag(D);
   zeroTol = 1e3*eps*max(abs(D));
   if any(D<-zeroTol)
      error(message('Control:analysis:getSectorIndex6'))
   end
   ip = find(D>zeroTol);
   M0 = Q;
   W1 = [];
   W2 = lrscale(U(:,ip),[],sqrt(D(ip)));
end

