function [F,LIST] = ffrf(THIS,FREQ,varargin)
% ffrf  Frequency response of transition variables to measurement variables.
%
% Syntax
% =======
%
%     [F,LIST] = ffrf(m,freq,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Model object for which the frequency response function
% will be computed.
%
% * `FREQ` [ numeric ] - Vector of frequencies for which the response
% function will be computed.
%
% Output arguments
% =================
%
% * `F` [ numeric ] - Array with frequency responses of transition
% variables (in rows) to measurement variables (in columns).
%
% * `LIST` [ cell ] - List of transition variables in rows of the `F`
% matrix, and list of measurement variables in columns of the `F` matrix.
%
% Options
% ========
%
% * `'maxIter='` [ numeric | *500* ] - Maximum number of iteration when
% computing the steady-state Kalman filter.
%
% * `'output='` [ *'namedmat'* | numeric ] - Output matrix `F` will be
% either a namedmat object or a plain numeric array; if the option
% `'select='` is used, `'output='` is always `'namedmat'`.
%
% * `'select='` [ char | cellstr | *`Inf`* ] - Return the frequency response
% function for selected variables only; `Inf` means all variables.
%
% * `'tolerance='` [ numeric | *1e-7* ] - Convergence tolerance when
% computing the steady-state Kalman filter.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

% Parse input arguments.
P = inputParser();
P.addRequired('m',@ismodel);
P.addRequired('freq',@isnumeric);
P.parse(THIS,FREQ);

% Parse options.
opt = passvalopt('model.ffrf',varargin{:});

if ischar(opt.select)
    opt.select = regexp(opt.select,'\w+','match');
end

if ischar(opt.exclude)
    opt.exclude = regexp(opt.exclude,'\w+','match');
end

isselect = iscellstr(opt.select);
isnamedmat = strcmpi(opt.output,'namedmat') || isselect;

% TODO: Implement the `'exclude='` option through the `'select='` option.

%**************************************************************************

ny = length(THIS.solutionid{1});
nx = length(THIS.solutionid{2});
nalt = size(THIS.Assign,3);

if ~isempty(opt.exclude)
    opt.exclude = regexprep(opt.exclude,'@?log\((.*?)\)','$1');
    exclude = strfun.findnames(THIS.name(THIS.nametype == 1),opt.exclude);
    index = isnan(exclude);
    if any(index)
        warning_(1,opt.exclude(index));
    end
    exclude(index) = [];
    aux = false([1,sum(THIS.nametype == 1)]);
    aux(exclude) = true;
    exclude = aux;
else
    exclude = false(1,ny);
end

FREQ = FREQ(:)';
nfreq = length(FREQ);
F = nan(nx,ny,nfreq,nalt);

if ny > 0
    [flag,nansol] = isnan(THIS,'solution');
    for ialt = find(~nansol)
        [T,R,k,Z,H,d,U,Omega] = mysspace(THIS,ialt,false);
        % Remove measurement variables excluded by user.
        Z(exclude,:) = [];
        H(exclude,:) = [];
        d(exclude,:) = [];
        % Compute FFRF.
        F(:,~exclude,:,ialt) = ...
            freqdom.ffrf2( ...
            T,R,k,Z,H,d,U,Omega, ...
            FREQ,opt.tolerance,opt.maxiter);
    end
    % Solution not available.
    if flag
        utils.warning('model', ...
            '#Solution_not_available',strfun.alt2str(nansol));
    end
end

% List of variables in rows and columns of `F`.
LIST = THIS.solutionvector(1:2);

% Convert output matrix to namedmat object.
if isnamedmat
    F = namedmat(F,LIST{2},LIST{1});
end

% Select requested variables.
if isselect
    F = select(F,opt.select);
end

end
