%
% Example 4.11, p.134 - Compute pole and zero directions
%
% G(s) = 1/(s+1) [ s-1    4    ]      (eq. 4.63)
%                [ 4.5  2(s-1) ]
%
g11 = nd2sys([1 -1], [1 2]);
g12 = nd2sys([1], [1 2],4);
g21 = nd2sys([1], [1 2],4.5);
g22 = nd2sys([1 -1], [1 2],2);

G = sysbal(sbs(abv(g11,g21),abv(g12,g22)))  % minimal realization has 1 state

---------------------------------------------------------------------
% First compute pole and zero directions the "crude" way using SVD
% This method is NOT generallay recommended, e.g., with several poles
% and zeros at the same location.
---------------------------------------------------------------------

p=spoles(G)        % one stable   (LHP) pole at p = -2
z=szeros(G)        % one unstable (RHP) zero at z = 4
pp = p*(1+1.e-4);  % move the pole a little away before evaluating G(p)
Gp = frsp(G,-1*pp);% Note that G(i*w)=frsp(w), so G(p) = G(i*-i*p)=frsp(-i*p)
[u,s,v]=vsvd(Gp);  % Take SVD of G(p)
yp = sel(u,':',1)  % Output pole direction is [0.5547  -0.8321]'
up = sel(v,':',1)  % Input  pole direction is [0.6     -0.8   ]'

Gg = vunpck(frsp(G,-i*4.76434)); r=rank(Gg) % r is Normal rank of G
Gz = frsp(G,-i*z); [u,s,v]=vsvd(Gz);
yz = sel(u,':',r)  % Output zero direction is [-0.8321 0.5547]'
uz = sel(v,':',r)  % Input  zero direction is [-0.8    0.6   ]'

---------------------------------------------------------------------
% Then compute the poles zeros from the state-space description.
% This is the PREFERRED way.
---------------------------------------------------------------------

% Pole directions first; see (4.69)
[A,B,C,D]=unpck(G); eig(A)   % the poles ...
[T,P] = eig(A);  YP = C*T    % associated pole output vectors
YPn = YP / norm(YP)          %   normalize to get pole output directions
eig(A')                      % the poles again ...
[Q,P] = eig(A'); UP = B'*Q   % associated pole input vectors
UPn= UP / norm(UP)           %   normalize to get pole input directions

% Now the zero directions, see (4.61) and (4.62)
% NOTE: Functions ozde and izde are available on the book's web-page
[Zy,YZ,Xy] = ozde(G)    % ozde: function for zero output directions
[Zu,UZ,Xu] = izde(G)    % izde: function for zero input directions


