function [INDX,w,FocusInfo,InfFlag] = sectorresp(D,~,W1,W2,wspec,isPlotted)
% Computes relative sector index as a function of frequency.
%
%  [INDX,W,FOCUSINFO] = sectorresp(D,M0,W1,W2,WSPEC,ISPLOTTED) computes
%  the sector index for a single MIMO model D over some user-defined or
%  auto-generated frequency grid W. For the R-index (M0=[]), INDX is of
%  size Nu-by-Nf and contains the generalized singular values of
%  (W1'*H(jw),W2'*H(jw)) at each frequency.
%
%  WSPEC specifies the frequency grid or range as follows:
%             [] :  none (auto-selected)
%    {fmin,fmax} :  user-defined frequency range (grid spans
%                   this range)
%         vector :  user-defined frequency grid
%
%  ISPLOTTED is true when the data is used for plotting and false
%  otherwise (interpolation is used only when ISPLOTTED=true).
%  See FREQFOCUS for details on the contents of FOCUSINFO.

%   Copyright 1986-2015 The MathWorks, Inc.
Ts = abs(D.Ts);
RealFlag = isreal(D);
[w,Focus] = parseFreqSpec(D,wspec,RealFlag && isPlotted);
AutoGrid = isempty(w);
AutoFocus = any(isnan(Focus));

% Scale the state-space realization to maximize accuracy in the
% frequency band of interest
if ~D.Scaled
    [D.a,D.b,D.c,D.e,~,~,Info] = ...
        ltipack.xscale(D.a,D.b,D.c,D.d,D.e,D.Ts,'Focus',Focus,'Warn',false);
    if ~isempty(Info.WarnID)
        % Tailor warning to SIGMA
        warning(message(strrep(Info.WarnID,...
            'Control:transformation:StateSpaceScaling',...
            'Control:analysis:ScalingIssue')))
    end
    D.Scaled = true;
end

% Check that H2(s) = W2'*H(s) is minimum phase (AFTER scaling)
InfFlag = ~hasSectorResp(D,W1,W2,Focus);
if InfFlag
   INDX = [];  w = [];  FocusInfo = [];  
   return
end

% Generate frequency grid if none specified
if AutoGrid
   % Get frequency grid from system dynamics (with all delays set to zero)
   if isequal(iosize(D),[1 1])
      % SISO
      [z,p] = iodynamics(D,'silent');
   else
      % MIMO: Use transmission zeros
      hw = ctrlMsgUtils.SuspendWarnings; %#ok<*NASGU>
      z = {tzero(D)};  p = {pole(D)};
      hw = [];
   end
   
   % Generate frequency grid
   OPTS = struct('Focus',Focus,'Negative',~RealFlag,...
      'Decades',AutoFocus && ~isPlotted);
   w = freqgrid(z,p,Ts,4,OPTS); % No delay-induced dynamics
else
   z = {zeros(0,1)};  p = {zeros(0,1)};  
end

% Compute indices
h = fresp(D,w);
INDX = ltipack.getINDEX(h,W1,W2);

% Pick frequency focus
FocusInfo = freqfocus(Focus,4,w,permute(h,[3 1 2]),z,p,Ts,0,[]);

% Sort frequencies if response is plotted
if isPlotted
    [w,is] = sort(w);  INDX = INDX(:,is);
end


