function iOPT = xscale_OptimGrid(w,g,relacc,iBP,slopes,wcf,Ts)
% Select grid for sensitivity minimization based on gain profile.
% Note: W always contains at least two frequencies.

%   Copyright 2020 The MathWorks, Inc.
nf = numel(w);
nBP = numel(iBP);
if nBP>0
   % Cluster break points W(iBP)
   GSEP = 1e3;  % min gain separation between clusters
   nc = 0;      % number of clusters
   is = zeros(nBP,1);  % start index
   ie = zeros(nBP,1);  % end index
   for ct=1:nBP
      j = iBP(ct);
      if nc==0 || g(j)>glast*GSEP || g(j)<glast/GSEP
         % Start new cluster
         nc = nc+1;  is(nc) = j;
      end
      ie(nc) = j;  glast = g(j);
   end
   
   % Pick focus cluster W(iOPT)
   if nc>1
      % A cluster is dominant if it has higher gain than its neighbors.
      % Find dominant clusters if any
      log2w = log2(w);   log2g = log2(g);
      cslopes = round([(log2g(is(1))-log2g(1))/(log2w(is(1))-log2w(1)) ; ...
         (log2g(ie(1:nc-1))-log2g(is(2:nc)))./(log2w(ie(1:nc-1))-log2w(is(2:nc))) ; ...
         (log2g(end)-log2g(ie(nc)))/(log2w(end)-log2w(ie(nc)))]);  % nc+1
      kdom = find(cslopes(1:nc)>0 & cslopes(2:nc+1)<0);
      ndom = numel(kdom);
      
      % Select "focus" cluster as:
      %   * Dominant cluster with highest gain if any cluster is dominant
      %   * Otherwise, score clusters to penalize eccentricity and steep
      %     slope, and pick cluster with lowest score.
      if ndom>0
         SCORE = zeros(ndom,1);
         for ct=1:ndom
            SCORE(ct) = mean(log2g(is(kdom(ct)):ie(kdom(ct)))); % average log-gain
         end
         [~,ibest] = max(SCORE);
         ibest = kdom(ibest);
      else
         SCORE = zeros(nc,1);
         % AS = average slope at each frequency point
         as = [slopes(1) ; (slopes(1:nf-2,:)+slopes(2:nf-1,:))/2 ; slopes(nf-1)];
         for ct=1:nc
            % Use cluster distance to WCF
            S1 = min(abs(log10(wcf)-log10(w([is(ct),ie(ct)]))));
            % Use mean slope to favor areas of flat gain
            S2 = abs(mean(as(is(ct):ie(ct))));
            %[S1 S2 S1+S2]
            SCORE(ct) = S1+S2;
         end
         [~,ibest] = min(SCORE);
      end
      iOPT = (is(ibest):ie(ibest)).';
   else
      iOPT = (is(1):ie(1)).';
   end
else
   % No breakpoints (slope is nearly constant). Focus on WCF
   wcf = max(w(1),min(wcf,w(nf)));
   if w(nf)>0.1*pi/Ts && abs(mean(slopes))>1.5
      % Discrete time with high slope: focus on Nyquist frequency
      % sys = ss(zpk([.1 .3i 0 0],[1 .99i 1 1 1 .6+.8i .6-.8i .3i .99i 0 0 0],1,1e12));
      iOPT = nf;  % w(nf) = safe Nyquist freq
   else
      % Use default center of focus and its neighbors (3 points)
      [~,iOPT] = min(abs(w-wcf));  % WCF could be a resonance missing from W
   end      
end
   
% Add points between DC and W(iOPT(1)) when gain flat to the left. Keep 
% an eye on relative accuracy: load dBoeingShuttle sys2
i0 = iOPT(1);
if all(abs(log10(g(1:i0-1,:)/g(i0)))<1) % gain flat
   while i0>1 && relacc(i0-1)<1
      i0 = i0-1;
   end
   iOPT = [(i0:iOPT(1)-1)';iOPT];
end

% Add point to the right when gain gently rolls off:
% load dNUM_Scaling A B C D; sys = dss(A,B,C,D,eye(size(A)));
i0 = iOPT(end);  MINACC = 1e3*relacc(i0);
if i0<nf && slopes(i0)>-1.5 && slopes(i0)<0 && relacc(i0+1)<MINACC
   iOPT = [iOPT;i0+1];
end

% Handle singleton (optimization requires at least two points)
if isscalar(iOPT)
   i0 = iOPT;
   if i0==1 || (i0<nf && g(i0+1)*relacc(i0+1)<g(i0-1)*relacc(i0-1))
      iOPT = [i0;i0+1];
   else
      iOPT = [i0-1;i0];
   end
end