%   [Po, Yp, Xpo, Spo] = opde(G,epp,RF)%%   The ``Output-Pole-Direction-Eigenvalue'' opde-function computes the %        output pole directions through the use of  egenvalue decomposition.%   NOTE:  Avoids the change in ordering of eigenvalues that may occur if %        the input and output directions are computed separately.%%   Outputs: Po   - Vector containing the poles.%            Yp   - Matrix containing the output pole vectors/directions (normalized = length 1), %                  column i corresponds to pole Po(i).%            Xpo   - Output state directions  = right eigenvectors of A (non-normalized).%                    NOTE: These are scaled such that pole vectors = pole directions%            Spo   - Scalars used to scale the state directions (eigenvectors).%                 %   Inputs:  G   - System matrix on mutools format: G=pck(A,B,C,D)%            epp - (optional) Tolerance, default value is 1e-12.%                  It is used to check the norm of the output directions.%                  If the norm of output direction is less than EPP,%                  zero is stored in the corresponding output direction.%            RF  - (optional) If FLG_RF==1 then the first element of %                  the output direction is made real.%%   Written by: Kjetil Havre, 12/9-1996, e-mail: kjetil@ife.no%   Requires subroutines: zerm.m and sorte.m%%   See also:  IPDE, PDFST, IZDE, OZDE, ZDSVD and PDSVD.%%   Reference: Havre K. and S. Skogestad, 1996, ``Effect of RHP Zeros and%              Poles on Performance in Multivariable Systems''.%   function [P,Y,X,S] = opde(G,epp,rf);   P = []; Y = []; X = [];    [mt,ny,nu, nx] = minfo(G);   S = ones(nx,1);      if strcmp(mt, 'syst')==0      disp( 'System matrix G is required, usage: [] = opde(G,epp)' );      return   end   if nargin < 3     rf = 0;   end   if nargin < 2      epp = 1e-12;   end   [A, B, C, D] = unpck(G);      if norm(D) > 1/epp      disp('Warning: Directions may be inaccurate due to large effect from D.');   end      [V, D] = eig(A,'nobalance');   P = diag(D);   [P,Is] = sorte(P); V = V(:,Is);   Hm = A*V-V*diag(P);   if norm(Hm) > 100*epp      disp('Warning: inaccurate eigenvalue computation.')      Hm = Hm   end      Npd = max(size(P));   for i=1:Npd	       X(:,i) = V(:,i);      if abs(P(i)-real(P(i))) < epp          P(i) = real(P(i));	 if isreal(X(:,i)) == 0 % Rotate eigenvector.	    inz = find(abs(X(:,i))>epp); 	    if isempty(inz) == 0     % Rotate eigenvector.	       vh = angle(X(inz(1),i));               X(:,i)=X(:,i)*exp(-sqrt(-1)*vh);	    end            X(:,i) = zerm(X(:,i),epp/10);         end      end            yph = C*X(:,i); nyph = norm(yph);            if nyph > epp         Y(:,i) = yph/nyph;	 X(:,i) = X(:,i)/nyph;	 	 S(i,1) = 1/nyph;	 if (isreal(Y(:,i))==0) & (rf==1)   % Rotate output direction.             inz = find(abs(Y(:,i))>epp); 	    if isempty(inz) == 0	       vh = angle(Y(inz(1),i));               Y(:,i) = Y(:,i)*exp(-sqrt(-1)*vh);                X(:,i)=X(:,i)*exp(-sqrt(-1)*vh);	       S(i,1) = S(i,1)*exp(-sqrt(-1)*vh);            end         end      else         Y(:,i) = zeros(ny,1);    % Norm is small unobservable.	 S(i,1) = 0;      end   end   P = zerm(P,epp/10); X = zerm(X,epp/10); Y = zerm(Y,epp/10);end