function [s,ndx]=dsort(p)
%DSORT  Sort complex discrete eigenvalues in descending order.
%
%   S = DSORT(P) sorts the complex eigenvalues in the vector P in 
%   descending order by magnitude.  The unstable eigenvalues (in
%   the continuous-time sense) will appear first.
%
%   [S,NDX] = DSORT(P) also returns the vector NDX containing the 
%   indexes used in the sort.
%
%   See also: ESORT and SORT.

%   Clay M. Thompson  7-23-90, AFP 6-1-94, PG 6-21-96
%   Copyright (c) 1986-97 by The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 1997/04/09 00:36:16 $

error(nargchk(1,1,nargin));
p = p(:);
if ~isequal(sort(p(imag(p)>0)),sort(conj(p(imag(p)<0)))),
   error('Complex eigenvalues must come in conjugate pairs.');
end

len = length(p);
s = zeros(len,1);

if nargout<=1,
   % Case NDX is not required: sort eigenvalues with nonnegative imaginary part
   p = p(imag(p)>=0);
   [trash,isort] = sort( -abs(p) );
   p = p(isort);

   k = 1;
   for e = p.',
      s(k) = e;   k = k+1;
      if imag(e),
         s(k) = conj(e);   k = k+1;
      end
   end
else
   % Same as above, but keep track of permutation NDX
   ndx = zeros(len,1);
   indp = find(imag(p)>=0);
   indn = find(imag(p)<0);
   [trash,isort] = sort( -abs(p(indp)) );
   indp = indp(isort);

   k = 1;
   for j = indp',
      e = p(j);
      s(k) = e;   ndx(k) = j;    k = k+1;
      if imag(e),
         jc = find(p(indn) == conj(e));
         jc = jc(1);
         s(k) = conj(e);   ndx(k) = indn(jc);    k = k+1;
         indn(jc) = [];
      end
   end
end

% end dsort
