home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / STRANS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.8 KB  |  70 lines
  1. % function [sysout,v] = strans(sys)
  2. %
  3. %   STRANS puts the A matrix in bidiagonal form
  4. %   with the complex conjugate roots in two by two form. 
  5. %   SYSOUT contains the transformed SYSTEM matrix and
  6. %   V is the transformation matrix.
  7. %
  8. %   See also: EIG, MMULT, REORDSYS, STATECC, SCLIN, 
  9. %             SCLOUT, and TRUNC.
  10.  
  11. function [sysout,vzi] = strans(sys);
  12.  if nargin ~= 1
  13.    disp('usage: [sysout,v] = strans(sys)')
  14.    return
  15.  end
  16.  [mtype,mrows,mcols,mnum] = minfo(sys);
  17.  if mtype == 'syst'
  18.    if mnum == 0
  19.      disp('input matrix has no states')
  20.    else
  21.     [a,b,c,d] = unpck(sys);
  22.     [v,dia] = eig(a);
  23. %  reorder by increasing magnitude
  24.     [dia,isort] = sort(diag(dia));
  25.     v = v(:,isort);
  26.     comp = imag(dia);
  27.     i = find(comp);
  28.     [n,m] = size(i);
  29.     nelem = max(n,m);
  30.     if ceil(nelem/2) ~= floor(nelem/2)
  31.       error(' there are an odd number of complex roots')
  32.       return
  33.     end
  34.     z = sqrt(1/2)*[1 sqrt(-1); sqrt(-1) 1];
  35.     old = 0;
  36.     mat = eye(a);
  37.     mat1 = mat;
  38. %    construct the transformation matrix
  39.     for ix = 1:nelem
  40.       if ix ~= old
  41.         mat(i(ix),i(ix)) = z(1,1);
  42.         mat(i(ix),i(ix)+1) = z(1,2);
  43.         mat(i(ix)+1,i(ix)) = z(2,1);
  44.         mat(i(ix)+1,i(ix)+1) = z(2,2);
  45.         mat1(i(ix),i(ix)) = 1;
  46.         mat1(i(ix),i(ix)+1) = 1;
  47.         mat1(i(ix)+1,i(ix)) = 1;
  48.         mat1(i(ix)+1,i(ix)+1) = 1;
  49.         old = ix+1;
  50.       else
  51.         old = 0;
  52.       end
  53.     end
  54.  
  55.     vzi = sign(real(v/mat)).*abs(v/mat);
  56.     vzii = sign(real(mat/v)).*abs(mat/v);
  57.     anew = vzii*a*vzi;
  58.     bnew = vzii*b;
  59.     cnew = c*vzi;
  60. %   zeros out some round off error
  61.     anew = anew.*mat1;
  62.     sysout = pck(anew,bnew,cnew,d);
  63.    end
  64.  else
  65.    error('input matrix is not a SYSTEM matrix')
  66.    return
  67.  end
  68. %
  69. % Copyright MUSYN INC 1991,  All Rights Reserved
  70.