function [alpha,w,FocusInfo] = diskmarginresp(L,sigma,wspec,STABLE,isPlotted)
% Generates disk margin data from open loop L and skew SIGMA.
%
%  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
%
%  STABLE indicates whether closed-loop model is stable. This is used
%  by TuningGoal.Margin to ensure consistency with SYSTUNE, EVALGOAL, etc.
%  When STABLE=[], stability is assessed directly from X.
%
%  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 2020 The MathWorks, Inc.
Ts = abs(L.Ts);
RealFlag = isreal(L);
[w,Focus] = parseFreqSpec(L,wspec,RealFlag && isPlotted);
AutoGrid = isempty(w);
AutoFocus = any(isnan(Focus));

% Skip this if already known to be closed-loop unstable
if ~isequal(STABLE,false)
   % Eliminate delays
   if hasdelay(L)
      if Ts==0
         L = pade(L,4,4,4);
      else
         L = elimDelay(L);
      end
   end
   
   % Form X = S+(sigma-1)/2*I = lft([(1+sigma)/2 -1;1 -1],L)
   [nL,~] = iosize(L);  % must be square
   M = kron([(1+sigma)/2 -1;1 -1],eye(nL));
   X = lft(createGain(L,M),L,nL+1:2*nL,nL+1:2*nL,1:nL,1:nL);
   
   % Scale the state-space realization to maximize accuracy in the
   % frequency band of interest
   if ~X.Scaled
      [X.a,X.b,X.c,X.e] = ...
         ltipack.xscale(X.a,X.b,X.c,X.d,X.e,X.Ts,'Focus',Focus,'Warn',false);
      X.Scaled = true;
   end
   
   % Compute dynamics of X and assess stability
   if AutoGrid
      if nL==1
         [z,p] = iodynamics(X,'silent');
      else
         % MIMO: Use transmission zeros
         hw = ctrlMsgUtils.SuspendWarnings; %#ok<*NASGU>
         z = {tzero(X)};  p = {pole(X)};
         hw = [];
      end
      if isempty(STABLE)
         STABLE = (Ts==0 && all(real(cat(1,p{:}))<0)) || ...
            (Ts>0 && all(abs(cat(1,p{:}))<1));
      end
   else
      z = [];  p = [];
      if isempty(STABLE)
         STABLE = isstable(X);
      end
   end
end


% Compute disk margin as a function of frequency
if STABLE
   % Generate frequency grid if none specified
   if AutoGrid
      % Get frequency grid from system dynamics (with all delays set to zero)
      OPTS = struct('Focus',Focus,'Negative',~RealFlag,...
         'Decades',AutoFocus && ~isPlotted);
      w = freqgrid(z,p,Ts,4,OPTS);
   end
   
   % Compute disk margin at grid points
   h = fresp(X,w);
   if nL==1
      UB = abs(h);
      if AutoGrid
         % Add peak
         [gpeak,fpeak] = norminf(X,1e-6);
         UB = [UB(:) ; gpeak];
         w = [w ; fpeak];
      end
   else
      bnds = mussv(h,ones(nL,2),'s');
      UB = bnds(1,1,:);
      if AutoGrid
         % Add peak
         [a,b,c,d] = getABCD(X);
         [blkData,opt] = ssmussvSetUp(ones(nL,2),'s',[]);
         [~,UBcert] = ssmussvPeak('robstab',a,b,c,d,Ts,blkData,opt);
         [UBmax,imax] = max([UBcert.gUB]);
         UB = [UB(:) ; UBmax];
         w = [w ; UBcert(imax).w];
      end
   end
   alpha = 1./UB(:);
   
   % Compute min margin
   if AutoGrid || isPlotted
      % Sort frequencies if response is plotted
      [w,is] = sort(w);  alpha = alpha(is,:);
   end
   
   % Pick frequency focus
   FocusInfo = freqfocus(Focus,4,w,alpha,z,p,Ts,0,[]);
else
   % Not stable: GM=1 and PM=0 at all frequencies
   if AutoGrid
      if AutoFocus
         w = logspace(-20,min(20,log10(pi/Ts)),10)';
      else
         w = logspace(log10(Focus(1)),log10(Focus(2)),10)';
      end
   end
   alpha = zeros(size(w));
   FocusInfo = struct('Focus',Focus,'DynRange',Focus,'Soft',AutoFocus);
end


