## Copyright (C) 1996, 2000, 2002, 2003, 2004, 2005, 2007
##               Auburn University.  All rights reserved.
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{a}, @var{b}, @var{c}, @var{d}] =} zp2ss (@var{zer}, @var{pol}, @var{k})
## Conversion from zero / pole to state space.
##
## @strong{Inputs}
## @table @var
## @item zer
## @itemx pol
## Vectors of (possibly) complex poles and zeros of a transfer
## function. Complex values must come in conjugate pairs
## (i.e., @math{x+jy} in @var{zer} means that @math{x-jy} is also in @var{zer}).
## The number of zeros must not exceed the number of poles.
## @item k
## Real scalar (leading coefficient).
## @end table
##
## @strong{Outputs}
## @table @var
## @item @var{a}
## @itemx @var{b}
## @itemx @var{c}
## @itemx @var{d}
## The state space system, in the form:
## @iftex
## @tex
## $$ \dot x = Ax + Bu $$
## $$ y = Cx + Du $$
## @end tex
## @end iftex
## @ifinfo
## @example
##      .
##      x = Ax + Bu
##      y = Cx + Du
## @end example
## @end ifinfo
## @end table
## @end deftypefn

## Author: David Clem
## Created: August 15, 1994

function [a, b, c, d] = zp2ss (zer, pol, k)

  if (nargin != 3)
    print_usage ();
  endif

  if (! (isvector (zer) || isempty (zer)))
    error ("zer(%d,%d) should be a vector", rows (zer), columns (zer));
  elseif (! (isvector (pol) || isempty (pol)))
    error ("pol(%d,%d) should be a vector", rows (pol), columns (pol));
  elseif (! isscalar(k))
    error ("k(%d,%d) should be a scalar", rows (k), columns (k));
  elseif (k != real (k))
    warning ("zp2ss: k is complex")
  endif

  zpsys = ss (zeros (0, 0), zeros (0, 1), zeros (1, 0), k);

  ## Find the number of zeros and the number of poles
  nzer = length (zer);
  npol = length (pol);

  if (nzer > npol)
    error ("%d zeros, exceeds number of poles=%d", nzer, npol);
  endif

  ## Sort to place complex conjugate pairs together
  zer = sortcom (zer);
  pol = sortcom (pol);

  ## construct the system as a series connection of poles and zeros
  ## problem: poles and zeros may come in conjugate pairs, and not
  ## matched up!

  ## approach: remove poles/zeros from the list as they are included in
  ## the ss system

  while (length (pol))

    ## search for complex poles, zeros
    cpol = [];
    czer = [];
    if (! isempty (pol))
      cpol = find (imag (pol) != 0);
    endif
    if (! isempty (zer))
      czer = find (imag (zer) != 0);
    endif

    if (isempty (cpol) && isempty (czer))
      pcnt = 1;
    else
      pcnt = 2;
    endif

    num = 1;      # assume no zeros left.
    switch (pcnt)
      case 1
	## real pole/zero combination
	if (length (zer))
	  num = [1, -zer(1)];
	  zer = zer(2:length(zer));
	endif
	den = [1, -pol(1)];
	pol = pol(2:length(pol));
      case 2
	## got a complex pole or zero, need two roots (if available)
	if (length (zer) > 1)
	  [num, zer] = __zp2ssg2__ (zer);       # get two zeros
	elseif (length (zer) == 1)
	  num = [1, -zer];                # use last zero (better be real!)
	  zer = [];
	endif
	[den, pol] = __zp2ssg2__ (pol);         # get two poles
      otherwise
	error ("pcnt = %d", pcnt);
      endswitch

      ## pack tf into system form and put in series with earlier realization
      zpsys1 = tf (num, den, 0, "u", "yy");

      ## change names to avoid warning messages from sysgroup
      zpsys  = syssetsignals (zpsys, "in", "u1", 1);
      zpsys1 = sysupdate (zpsys1, "ss");
      nn     = sysdimensions (zpsys);        # working with continuous system
      zpsys  = syssetsignals (zpsys, "st", __sysdefioname__ (nn, "x"));
      nn1    = sysdimensions (zpsys1);
      zpsys1 = syssetsignals (zpsys1, "st", __sysdefioname__ (nn1, "xx"));

      zpsys = sysmult (zpsys, zpsys1);

    endwhile

    [a, b, c, d] = sys2ss (zpsys);

endfunction