function [f, g] = sigmin(x, pars)
% call:  [f, g] = sigmin(x, pars)
% function variables are x(1), x(2) (real and imag parts of a complex no z)
% pars.A and pars.B are specified matrices
% corrected October 2017 because there was a bug for complex matrices
% previous definition:  sigmin(x; A,B) = sigma_min([zI - A' ])
%                                                  [   B'   ]
% but this does convenient transposition does not work if A is complex. 
% Using this transposition, which has its conveniences, requires also
% conjugating z:
%                       sigmin(x; A,B) = sigma_min([conj(z)I - A' ])
%                                                  [         B'   ] 

% The minimum value of this function over all z is
% the distance from the control system
%
%                dy/dt = Ay + Bu
%
% to the nearest uncontrollable system (a system which cannot be controlled
% from any given starting state y_0 to any given finishing state y_finish).
% Here A is n by n, B is n by p, y is an n-vector of states and u is
% a p-vector of control inputs.  
% References: Trefethen and Wright, IMA J Num Anal, 2002 
%             Eising, Sys Control Letters, 1984
% Added in 2017: did not know about Paige, 1981, Trans Auto Control
% which learned about from nice Master's thesis by Boyce at Va. Tech, 2010
%
n = size(pars.A, 1);
if n ~= size(pars.A, 2)
   error('pars.A must be square')
end
if n ~= size(pars.B, 1)
   error('pars.A and pars.B must have same row dimension')
end
p = size(pars.B, 2);
if p >= n
   error('pars.B cannot have more columns than rows')
end
z = x(1) + 1i*x(2);
M = [conj(z)*eye(n) - (pars.A)';  (pars.B)']; % previously no conj, corrected 2017
[U, S, V] = svd(M,0);   % economy size
f = S(n,n);             % smallest nontrivial singular value
u = U(:,n);             % has length n+p
v = V(:,n);             % has length p
uchopped = u(1:n);      % to give us u'[ I ] v
w = uchopped'*v;        %              [ 0 ] 
g = [real(w); imag(w)]; % previously -imag(w), corrected 2017
