% function out = frsp(sys,omega,T,balflg)
%
%   Complex frequency response of a SYSTEM matrix, produces
%   a VARYING matrix which contains its frequency response:
%
%    SYS     - SYSTEM matrix
%    OMEGA   - response calculated at these frequencies, can
%                also be another VARYING matrix, and then these
%                independent variables are used
%    T       - 0 (default) continuous system, if T ~= 0 then
%		  discrete evaluation with sample time T is used.
%    BALFLG  - 0 (default) balances A matrix prior to eval, nonzero
%                 value for BALFLG leaves state-space unchanged
%    OUT     - VARYING frequency response matrix
%
%   See also: BALANCE, HESS, SAMHLD, TUSTIN, and VPLOT.

function out = frsp(sys,omega,T,balflg)
  if nargin < 2
    disp('usage: out = frsp(sys,omega)')
    return
  elseif nargin == 2
   balflg = 0;
   T = 0;
  elseif nargin == 3
   balflg = 0;
  end
  [mtype,mrows,mcols,mnum] = minfo(sys);
  if mtype == 'vary'
    error('input matrix is already a VARYING matrix')
    return
  elseif mtype == 'empt'
    out = [];
    return
  end
  [omtype,omrows,omcols,omnum] = minfo(omega);
  if omtype == 'cons'
    if omrows == 1
      omega = omega.';
    end
    npts = length(omega);
  elseif omtype == 'vary'
    omega = omega(1:omnum,omcols+1);
    npts = omnum;
  else
    error('omega should be a vector, or VARYING matrix')
    return
  end
  if mtype == 'cons'
    out == [];
    for i=1:npts
      out = [out; sys];
    end
    out = [out [omega;zeros(npts*(mrows-1),1)];zeros(1,mcols-1) npts inf];
    return
  end
  if isempty(omega)
    out = [];
    return
  else
    comega = sqrt(-1) * omega;
    if T ~= 0
      comega = exp(T*comega);
    end
    [a,b,c,d] = unpck(sys);
    if nargin == 2 | balflg == 0
      [t,a] = balance(a);
      b = t\b;
      c = c*t;
    end
    [states dum]=size(a);
    [outputs inputs] = size(d);
    [p,hesa] = hess(a);
    hesb = p' * b;
    hesc = c * p;
    npts = length(comega);
    idn = eye(states);
    nrout = mrows;
    ncout = mcols;
    out = zeros((nrout*npts)+1,ncout+1);
    fout = (npts+1)*nrout;
    pout = 1:nrout:fout;
    poutm1 = pout(2:npts+1) - 1; 
    for i=1:npts
      g = (comega(i) * idn - hesa) \ hesb;
      out(pout(i):poutm1(i),1:ncout) = d + hesc*g;
    end
    out(1:npts,ncout+1) = omega;
    out((nrout*npts)+1,ncout+1) = inf;
    out((nrout*npts)+1,ncout) = npts;
  end
%
% Copyright MUSYN INC 1991,  All Rights Reserved
