function [resp,freq,Ts,covresp] = frdata(sys,varargin)
%FRDATA  Quick access to frequency response data.
%
%   [RESP,FREQ] = FRDATA(SYS) returns the response data RESP and frequency 
%   vector FREQ of the FRD model SYS. The frequencies FREQ are expressed 
%   in the units specified in SYS.FrequencyUnit.
%
%   For a single model SYS with Ny outputs, Nu inputs, and Nw frequency
%   points, RESP is an array of size [Ny Nu Nw] where the (I,J,K) element 
%   specifies the response from input J to output I at the K-th frequency 
%   point FREQ(K).  
%
%   [RESP,FREQ,TS] = FRDATA(SYS) also returns the sample time TS. Other
%   properties of SYS can be accessed using struct-like dot syntax (for
%   example, SYS.FrequencyUnit).
%
%   [RESP,FREQ,TS,COVRESP] = FRDATA(SYS) also returns the covariance
%   COVRESP of the response data RESP for IDFRD model SYS. COVRESP is a
%   5D-array where covH(i,j,k,:,:) contains the 2-by-2 covariance matrix of
%   the response RESP(i,j,k).  The (1,1) element is the variance of the
%   real part, the (2,2) element the variance of the imaginary part and the
%   (1,2) and (2,1) elements the covariance between the real and imaginary
%   parts.
%
%   For SISO models, the syntax
%       [RESP,FREQ,TS,COVRESP] = FRDATA(SYS,'v')
%   returns RESP as a column vector rather than a 3-dimensional array. For
%   IDFRD model SYS, COVRESP is returned as a 3D, rather than 5D, array of
%   size NF-by-2-2, where NF = length(FREQ).
%
%   [RESP,FREQ,TS] = FRDATA(SYS,J1,...,JN) extracts the data for the 
%   (J1,...,JN) entry in the model array SYS.
%
%   See also FRD, GENFRD, IDFRD.

%   Author(s): P. Gahinet, S. Almy
%   Copyright 1986-2011 The MathWorks, Inc.

% Detect 'v' flag
% String support: strncmpi does not support cells for strings
vflag = find(cellfun(@(x) strncmpi(x,'v',1),varargin));
varargin(:,vflag) = [];

% Get data
try
   if isempty(varargin)
      % Possibly multiple models
      s = size(sys);
      nsys = prod(s(3:end));
      if nsys==0
         resp = zeros([s(1:2) 0 s(3:end)]);
         freq = zeros(0,1);
         Ts = 0;
         covresp = [];
      else
         nf = size(sys,'freq');
         resp = zeros([s(1:2) nf s(3:end)]);
         nu = s(2);
         NoCov = true;
         if isa(sys,'idlti')
            if nu==0
               resp = zeros([s(1) s(1) nf s(3:end)]);
               covresp = NaN([s(1) s(1) nf s(3:end)]);
            else
               covresp = NaN([s(1:2) nf 2 2 s(3:end)]);
            end            
         else
            covresp = [];
         end
         
         for ct = 1:nsys
            [resp(:,:,:,ct),freq,Ts,covresp_] = frdata_(sys,ct);
            if ~isempty(covresp_)
               NoCov = false;
               if nu==0
                  covresp(:,:,:,ct) = covresp_;
               else
                  covresp(:,:,:,:,:,ct) = covresp_;
               end
            end
         end
         
         if NoCov, covresp = []; end            
      end
   else
      % Single model
      nsys = 1;
      [resp,freq,Ts,covresp] = frdata_(sys,varargin{:});
   end
catch ME
   ltipack.throw(ME,'command','frdata',class(sys))
end

% Return vectors if requested
if ~isempty(vflag) && nsys==1 && issiso(sys)
   % Convenience syntax for SISO case
   resp = resp(:);
   if ~isequal(covresp,[])
      covresp = permute(covresp,[3 4 5 1 2]);
   end
end
