function [WCG,Info,IRD] = wcgain(D,Opt)
% [StabMarg,Info,IRData] = WCGAIN(D,ULEVEL,OPT)
%     D      Scalar lftdataM object
%     OPT    rctoptions.wcAnalysis object

%   Copyright 2010-2015 The MathWorks, Inc.
[M,B,~,muB] = ulftdata(D);
if ~muB.ismuStandard
   % Structures with repeated, non-scalar blocks are not supported
   error(message('Robust:analysis:OneIllegalRSRP'))  % message OK for gain too
end
M = double(M); % M has no blocks, but is still LFT
BlockNames = {muB.BlkInfo.BlockName}';
nblk = numel(BlockNames);
ulevel = Opt.ULevel;

if nblk==0
   g = norm(M);
   WCG = struct('LowerBound',g,...
      'UpperBound',g);
   WCU = struct;
   Info = struct(...
      'Model',[],...
      'WorstPerturbation',WCU,...
      'Sensitivity',struct);
else
   usqrt = sqrt(ulevel);
   [ny,nu] = iosize(B);
   M(:,1:ny,:) = M(:,1:ny,:)*usqrt;
   M(1:nu,:,:) = M(1:nu,:,:)*usqrt;
   [ne,nd] = iosize(D);
   wcgainBlk = [muB.muBlkNby2;nd ne];
   fixedBlkIdx = 1:size(muB.muBlkNby2,1);  % uncertain blocks are "fixed" in WCG
   [bnds,tinfo] = mussv(M,wcgainBlk,Opt.MussvOptions,fixedBlkIdx);

   % Populate WCG
   WCG = struct(...
      'LowerBound',bnds(2),...
      'UpperBound',bnds(1));
   
   % Populate INFO
   PertM = mussvunwrap(tinfo.pvec,tinfo.blk);
   PertM = PertM(1:ny,1:nu); % keep only actual perturbations
   WCU = getWorstCasePerturbation(B,muB,ulevel*PertM,NaN);
   Info = struct(...
      'Model',[],...
      'WorstPerturbation',WCU,...
      'Sensitivity',[]);
   
   % Compute sensitivity data
   if strcmp(Opt.Sensitivity,'on')
      SensValues = localComputeSensitivity(M,wcgainBlk,bnds,Opt);
   else
      SensValues = NaN(nblk,1);
   end
   Info.Sensitivity = cell2struct(num2cell(SensValues),BlockNames,1);
end

% Report data
IRD = [];


%---------- Local subroutine ------------------------------%


function SensValues = localComputeSensitivity(M,wcgainBlk,bnds,Opt)
% Computes sensitivity of MUSSV peak
gamma = Opt.SensitivityPercent/100;
nblk = size(wcgainBlk,1);
blkData = rctutil.mkBlkData(wcgainBlk);
sfactor = zeros(nblk-1,1);
pku0 = bnds(1,1);
% Options for Sensitivity should be Silent, and take lowerBound down
% to only 1 try (so multipleTimes = 0).
SensMuOpt = [Opt.MussvOptions 's'];
mloc = find(SensMuOpt=='m');
if ~isempty(mloc)
   SensMuOpt(mloc+1) = '0';
end
for k = 1:nblk-1  % don't compute sensitivity wrt performance block
   tmpmat = M;
   affectedRows = blkData.ridxm(k):blkData.ridxm(k+1)-1;
   tmpmat(affectedRows,:,:) = (1+gamma)*tmpmat(affectedRows,:,:);
   sbnds = mussv(tmpmat,wcgainBlk,SensMuOpt,1:nblk-1);
   pku = sbnds(1,1); % 1-by-1-by-(1,2 or more)
   if pku0==0 && pku>0
      sfactor(k) = inf;
   elseif pku0==0 && pku==0
      sfactor(k) = 1;
   else
      sfactor(k) = abs(pku-pku0) / (gamma*pku0);
   end
end
SensValues = round(100*sfactor);
