## Copyright (C) 1996, 1998, 1999, 2000, 2002, 2003, 2004, 2005, 2006,
##               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{sys} =} sysgroup (@var{asys}, @var{bsys})
## Combines two systems into a single system.
##
## @strong{Inputs}
## @table @var
## @item asys
## @itemx bsys
## System data structures.
## @end table
##
## @strong{Output}
## @table @var
## @item sys
## @math{sys = @r{block diag}(asys,bsys)}
## @end table
## @example
## @group
##          __________________
##          |    ________    |
## u1 ----->|--> | asys |--->|----> y1
##          |    --------    |
##          |    ________    |
## u2 ----->|--> | bsys |--->|----> y2
##          |    --------    |
##          ------------------
##               Ksys
## @end group
## @end example
## The function also rearranges the internal state-space realization of @var{sys}
## so that the continuous states come first and the discrete states come last.
## If there are duplicate names, the second name has a unique suffix appended
## on to the end of the name.
## @end deftypefn

## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
## Created: August 1995
## modified by John Ingram July 1996
## A. S. Hodel: modified for variable number of arguments 1999

function sys = sysgroup (varargin)

  if (nargin < 1)
    print_usage ();
  endif

  ## collect all arguments
  arglist = {};
  for kk = 1:nargin
    arglist(kk) = varargin{kk};
    if (! isstruct (arglist{kk}))
      error ("sysgroup: argument %d is not a data structure", kk);
    endif
  endfor

  if (nargin == 2)
    ## the usual case; group the two systems together
    Asys = arglist{1};
    Bsys = arglist{2};

    ## extract information from Asys, Bsys to consruct sys
    Asys = sysupdate (Asys, "ss");
    Bsys = sysupdate (Bsys, "ss");
    [n1, nz1, m1, p1] = sysdimensions (Asys);
    [n2, nz2, m2, p2] = sysdimensions (Bsys);
    [Aa, Ab, Ac, Ad, Atsam, An, Anz, Ast, Ain, Aout, Ayd] = sys2ss (Asys);
    [Ba, Bb, Bc, Bd, Btsam, Bn, Bnz, Bst, Bin, Bout, Byd] = sys2ss (Bsys);
    nA = An + Anz;
    nB = Bn + Bnz;

    if (p1*m1*p2*m2 == 0)
      error ("sysgroup: argument lacks inputs and/or outputs");

    elseif (Atsam + Btsam > 0 && Atsam * Btsam == 0)
      warning ("sysgroup: creating combination of continuous and discrete systems")

    elseif (Atsam != Btsam)
      error ("sysgroup: Asys.tsam=%e, Bsys.tsam =%e", Atsam, Btsam);
    endif

    if (nA*nB > 0)
      A12 = zeros (nA, nB);
    else
      A12 = [];
    endif
    A = [Aa, A12; A12', Ba];

    if (nA*m2 > 0)
      B12 = zeros (nA, m2);
    else
      B12 = [];
    endif
    if (nB*m1 > 0)
      B21 = zeros (nB, m1);
    else
      B21 = [];
    endif
    if (isempty(Ab))
      Ab = [];
    endif
    if (isempty (Bb))
      Bb = [];
    endif
    B = [Ab, B12; B21, Bb];

    if (p1*nB > 0)
      C12 = zeros (p1, nB);
    else
      C12 = [];
    endif
    if (p2*nA > 0)
      C21 = zeros (p2, nA);
    else
      C21 = [];
    endif
    C = [Ac, C12; C21, Bc];

    if (p1*m2 > 0)
      D12 = zeros (p1, m2);
    else
      D12 = [];
    endif
    if (p2*m1 > 0)
      D21 = zeros (p2, m1);
    else
      D21 = [];
    endif
    D = [Ad, D12; D21, Bd];
    tsam = max (Atsam, Btsam);

    ## construct combined signal names; stnames must check for pure gain blocks
    if (isempty (Ast))
      stname = Bst;
    elseif (isempty (Bst))
      stname = Ast;
    else
      stname= __sysconcat__ (Ast, Bst);
    endif
    inname = __sysconcat__ (Ain, Bin);
    outname = __sysconcat__ (Aout, Bout);

    ## Sort states into continuous first, then discrete
    dstates = ones (1, (nA+nB));
    if (An)
      dstates(1:(An)) = zeros (1, An);
    endif
    if (Bn)
      dstates((nA+1):(nA+Bn)) = zeros (1, Bn);
    endif
    [tmp, pv] = sort (dstates);
    A = A(pv,pv);
    B = B(pv,:);
    C = C(:,pv);
    stname = stname (pv);

    ## check for duplicate signal names
    inname = __sysgroupn__ (inname, "input");
    stname = __sysgroupn__ (stname, "state");
    outname = __sysgroupn__ (outname, "output");

    ## mark discrete outputs
    outlist = find ([Ayd, Byd]);

    ## build new system
    sys = ss (A, B, C, D, tsam, An+Bn, Anz+Bnz, stname, inname, outname);

  else
    ## multiple systems (or a single system); combine together one by one
    sys = arglist{1};
    for kk = 2:length(arglist)
      printf ("sysgroup: kk=%d\n", kk);
      sys = sysgroup (sys, arglist{kk});
    endfor
  endif

endfunction