function [cl,w] = cloci(Z1,Z2,Z3,Z4,Z5,Z6)
%CLOCI  Characteristic Loci frequency response of continuous systems.
%	CLOCI(A,B,C,D) or CLOCI(SS_) produces a Characteristic Loci
%       plot of the complex matrix      
%                                            -1
%                             G(jw) = C(jwI-A) B + D 
%
%	The Char. loci are an extension of the Nyquist response for MIMO 
%	systems.  The Char. loci are given by the eigenvalues lambda_i(G(jw))
%	for 1<=i<=n where n is the size of G(jw), as jw traverses the Nyquist
%	D contour. Typically G(jw) will represent the interconnection of a 
%	controller and a plant. The frequency range and number of points are 
%	chosen automatically.  For square systems, CLOCI(A,B,C,D,'inv') produces 
%	the Char. Gain/Phase of the inverse complex matrix
%	     -1               -1      -1
%           G (jw) = [ C(jwI-A) B + D ]
%
%	CLOCI(A,B,C,D,W) or CLOCI(A,B,C,D,W,'inv') uses the user-supplied
%	frequency vector W which must contain the frequencies, in 
%	radians/sec, at which the Char. loci is to be 
%	evaluated. When invoked with left hand arguments,
%		[CL,W] = CLOCI(A,B,C,D,...)
%	returns the frequency vector W and the matrix CL with as many 
%	columns	as MIN(NU,NY) and length(W) rows, where NU is the number
%	of inputs and NY is the number of outputs.  No plot is drawn on 
%	the screen.  The Char. Loci are returned in descending order.
%
%       NOTICE: Modified by STC to return the complex-valued loci itselt rather
%       than the absolute value and phase of the loci. The loci are plotted as
%       real and imaginary values in a Nyquist-curve-like fashion.
%
%	See also: CGLOCI,LOGSPACE,SEMILOGX,NICHOLS,NYQUIST and BODE.

% S. Toffner-Clausen. Revised:  95.11.01, 99.31.03 (PA)
% Copyright (c) by the author
% All Rights Reserved.
nag1=nargin;
if nargin==0, eval('exresp(''cloci'',1)'), return, end

inargs='(a,b,c,d,w,invflag)';
eval(mkargs(inargs,nargin,'ss'))

if nag1==6,      % Trap call to RCT function
  if ~isstr(invflag)
    eval('[cg,ph] = cgloci2(a,b,c,d,w,invflag);')
    return
  end
  if ~length(invflag)
    nag1 = nag1 - 1
  end
end
nargin

error(nargchk(4,6,nag1));
error(abcdchk(a,b,c,d));

% Detect null systems
if ~(length(d) | (length(b) &  length(c)))
       return;
end

% Determine status of invflag
if nag1==4, 
  invflag = [];
  w = [];
elseif (nag1==5)
  if (isstr(w)),
    invflag = w;
    w = [];
    [ny,nu] = size(d);
    if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
  else
    invflag = [];
  end
else
  [ny,nu] = size(d);
    if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
end

% Generate frequency range if one is not specified.

% If frequency vector supplied then use Auto-selection algorithm
% Fifth argument determines precision of the plot.
if ~length(w)
  w=freqint(a,b,c,d,30);
end

[nx,na] = size(a);
[no,ns] = size(c);
nw = max(size(w));

% Balance A
[t,a] = balance(a);
b = t \ b;
c = c * t;

% Reduce A to Hessenberg form:
[p,a] = hess(a);

% Apply similarity transformations from Hessenberg
% reduction to B and C:
b = p' * b;
c = c * p;

s = w * sqrt(-1);
I=eye(length(a));
if nx > 0,
  for i=1:length(s)
    if isempty(invflag),
      char = c*((s(i)*I-a)\b) + d;
      cl(:,i) = eig(char);
    else
      char  = inv(c*((s(i)*I-a)\b) + d);
      cl(:,i) = eig(char);
    end
  end
else
  for i=1:length(s)
    if isempty(invflag)
      cl(:,i) = eig(d);
    else
      cl(:,i) = eig(inv(d));
    end
  end
end
  
% If no left hand arguments then plot graph.
if nargout==0
  subplot(111)
  plot(real(cl)',imag(cl)','-');
  ylabel('Imag. axis')
  xlabel('Real Axis');
  return % Suppress output
end