function wc = getGainCrossover(D,g)
% Crossover frequencies for specific gain level.

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

% Handle internal delays
if norm_hasInternalDelay(D)
   % Internal delays -> use interpolation
   [mag,phase,w] = freqresp(D,4,[],false);
   h = permute(mag.*exp(complex(0,phase)),[2 3 1]);
   Dfrd = ltipack.frddata(h,w,Ts);
   wc = getGainCrossover(Dfrd,g);
   return
end

% Simplify realization
D = sminreal(D);

% Compute gain crossings
[a,b,c,d,~,e] = getABCDE(D);
if isempty(a)
   % No states
   wc = zeros(0,1);  return
end

% Scale data
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);
   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);
   end
end

% Perform Hamiltonian test
% Note: Eigenvalues with magnitude less than wZERO may be inaccurate
[heigs,w_acc] = HamGainTest(a,b,c,d,e,Ts,g);

% 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
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)
   wc = wc(wc>=0);
end

% Evaluate gain at candidate frequencies
% Note: Some crossovers may be associated with non-maximal singular
% values so validation is always necessary
if Ts>0
   s = exp(complex(0,Ts*wc));
else
   s = complex(0,wc);
end
hws = ssfresp(a,b,c,d,e,s);
nf = numel(wc);
isOK = true(nf,1);
if isinf(g)
   for ct=1:nf
      isOK(ct) = isinf(norm(hws(:,:,ct)));
   end
else
   for ct=1:nf
      isOK(ct) = abs(1-norm(hws(:,:,ct))/g)<0.01;
   end
end
wc = unique(wc(isOK,:));