## -*-octave-*-  
##
## k = dcgain(sys)
##
## Compute the DC gain k of a linear time invariant (LTI) system
##
## Transfer function for a continuous state space system (A,B,C,D)
## G(s) = C * inv(s*I - A) * B + D
##
## DC Gain: evaluate G(s) as s -> 0
## k = C * inv(-A) * B + D
##
## Transfer function for a discrete state space system (A,B,C,D,T)
## G(z) = C * inv(z*I - A) * B + D
##
## DC Gain: evaluate G(z) as z -> 1
## k = C * inv(I-A) * B + D

## Copyright (C) 2016 Geraint Paul Bevan <geraint.bevan@gcu.ac.uk>
##
## This program 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.
##
## This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

## -*-texinfo-*-
## @deftypefn {Function File} {@var{sys} =} dcgain (@var{sys})
## Compute the dcgain of a linear time invariant (LTI) system.
##
## @strong{Inputs}
## @table @var
## @item sys
## A LTI model created by tf(), ss(), dss(), etc.
## @end table
##
## @strong{Outputs}
## @table @var
## @item k
## DC Gain of the LTI model
## @end table
##
## @strong{Example}
## @example
## G = Transfer function 'G' from input 'u1' to output ...
##
##                1          
## y1:  ---------------------
##      s^3 + 2 s^2 + 3 s + 4
##
##
## octave:1> k = dcgain(G)
##
## k =  0.25000
## @end example
##
## @seealso{tf,ss,dss}

function k = dcgain(sys)

  if (! isa(sys,"lti"))
    error("dcgain: argument must be an lti system");
  endif
  
  tsam = get(sys,'tsam');
  
  if (0 == tsam)		# Continuous system
    k = sys.c*inv(-sys.a)*sys.b+sys.d
  else				# Discrete system
    k = sys.c*inv(eye(size(sys.a))-sys.a)*sys.b+sys.d
  end

endfunction
