function wc = getSectorCrossover(D,W1,W2)
% Crossover frequencies for sector bound.

%   Copyright 1986-2011 The MathWorks, Inc.
Ts = abs(D.Ts);
JWTOL = 10*sqrt(eps);

% Can't handle delays in continuous time
if hasdelay(D)
   if Ts==0
      % Use interpolation
      [mag,phase,w] = freqresp(D,4,[],false);
      h = permute(mag.*exp(complex(0,phase)),[2 3 1]);
      wc = getSectorCrossover(ltipack.frddata(h,w,Ts),W1,W2);
      return
   else
      D = elimDelay(D,D.Delay.Input,D.Delay.Output,D.Delay.Internal);
   end
end

% Compute gain crossings
D = sminreal(D);
[a,b,c,d,~,e] = getABCDE(D);
if isempty(a)
   % No states
   wc = zeros(0,1);  return
end
if ~D.Scaled
   % Rescale for freq. response accuracy
   [a,b,c,e] = ltipack.xscale(a,b,c,d,e,Ts,'Warn',false);
end

% In CT, eliminate cancellations at Inf because frequency response can't be 
% reliably evaluated near s=Inf when E is singular
if Ts==0 && ~isempty(e)
   Ds = ltipack.ssdata(a,b,c,d,e,0);  % D may contain residual internal delays
   Ds.Scaled = true;
   [isP,Ds] = isproper(Ds);
   if isP
      % System is proper. Use reduced model with nonsingular E
      [a,b,c,d,~,e] = getABCDE(Ds);
   else
      error(message('Control:analysis:getSectorIndex10'))
   end
end

% Check G has normal rank
[ny,nu] = size(d);
if ~isequal(d(ny-nu+1:ny,:),eye(nu)) && normrank(a,b,c,d,e,1e3*eps)<nu
   error(message('Control:analysis:getSectorIndex11'))
end

% Perform Hamiltonian test for R=1
% Note: Eigenvalues with magnitude outside W_ACC may be inaccurate
[heigs,w_acc] = HamConicTest(a,b,c,d,e,Ts,[],W1,W2,-1,1);

% Pure imaginary eigenvalues provide candidate crossover frequencies
wn = abs(heigs);
wc = imag(heigs(abs(real(heigs))<JWTOL*wn | wn<=w_acc(1) | wn>=w_acc(2),:));

% Add WINF (excluded from Hamiltonian test)
% Note: Beware of introducing duplicate crossing near wInf
wInf = pi/Ts;
wc = [-wInf ; wc(abs(wc)<.9999*wInf) ; wInf];

% Keep only nonnegative frequencies for real-valued systems
if isreal(a) && isreal(b) && isreal(c) && isreal(d) && ...
      isreal(e) && isreal(W1) && isreal(W2)
   wc = wc(wc>=0);
end

% Evaluate R-index at candidate frequencies
% Note: Some crossovers may be associated with non-maximal R values 
% or with zeros of H(jw) so validation is always necessary
if Ts>0
   s = exp(complex(0,Ts*wc));
else
   s = complex(0,wc);
end
W12 = [W1 W2];
[h,~,~,hscale] = ssfresp(a,b,W12'*c,W12'*d,e,s);
nf = numel(wc);
nw1 = size(W1,2);
nw = size(W12,2);
isOK = true(nf,1);
for ct=1:nf
   R = ltipack.util.gsvmax(h(1:nw1,:,ct),h(nw1+1:nw,:,ct),hscale(ct));
   isOK(ct) = abs(1-R)<0.01;
end
wc = unique(wc(isOK,:));