function [stabmarg,DSV,info,SensText] = rsEngine(M,B,muB,Opt,Need)
% Inputs:
%   * M is lftdataFRD without blocks, so M.IC is frddata
%   * B is LFTBlockWrapper
%   * MUB is a description of the block structure
%   * OPT is a rctoptions.robuststab object
%   * NEED is a STRUCT with 4 fields: Report, Mussv4Info,
%     AllBadUncertainValues, Sensitivity
%
% Outputs:
%   * STABMARG: STRUCT with scalar-valued fields LowerBound, UpperBound,
%     DestabilizingFrequency
%   * DSV: STRUCT (destabilizing value for each block)
%   * INFO: STRUCT containing the following data:
%       INFO.Sensitivity: STRUCT (fields = block names)
%       INFO.Frequency: Same as M.IC.Frequency
%       INFO.BadUncertainValues: NF-by-1 struct when Need.AllBUVs,
%          empty otherwise (NF = number of frequency points)
%       INFO.MussvBnds: 1-by-2 DOUBLE
%       INFO.MussvInfo: 1-by-1 STRUCT (empty if ~Need.Mussv4Info)
%   * SENSTEXT is CHAR array

%   Author(s): MUSYN
%   Copyright 2010-2011 The MathWorks, Inc.

% Default outputs
stabmarg = struct('LowerBound',NaN,'UpperBound',NaN,...
   'DestabilizingFrequency',NaN);
info = struct('Sensitivity',struct,...
   'Frequency',NaN,...
   'BadUncertainValues',{cell(0,1)},...
   'MussvBnds',NaN(1,2),...
   'MussvInfo',struct);

if isempty(muB.BlkInfo)
   % No blocks
   [stabmarg,DSV,info,SensText] = rctutil.fillNoBlocks(Need,M,B,muB);      
else
   % Grab TS and Frequency from M.IC
   BlockName = {muB.BlkInfo.BlockName}';
   Frequency = M.IC.Frequency;
   nblk = numel(muB.BlkInfo);
   npts = length(Frequency);
   info.Frequency = Frequency;
   
   % Absorb all Delays, but do not interpolate.  It is true that FRD-like
   % objects never have internal delays.   It is also true that the InputDelay
   % and OutputDelay will always be identically zero on "Block" channels.
   % Nevertheless, FRESP is the simplest manner in which to obtain the 3-d
   % DOUBLE needed for the MUSSV computation.
   Mdata = fresp(M.IC,Frequency);
   szB = iosize(B);
   M11 = Mdata(1:szB(2),1:szB(1),:);  % ULFTDATA [Blocks;IO], M11, 3-D double-array
   [bnds,tinfo] = mussv(M11,muB.muBlkNby2,Opt.Mussv);
   if Need.Mussv4Info
      Ts = M.IC.Ts;
      MussvInfo = tinfo;
      MussvInfo.bnds = ltipack.frddata(MussvInfo.bnds,Frequency,Ts);
      MussvInfo.dvec = ltipack.frddata(MussvInfo.dvec,Frequency,Ts);
      MussvInfo.pvec = ltipack.frddata(MussvInfo.pvec,Frequency,Ts);
      MussvInfo.gvec = ltipack.frddata(MussvInfo.gvec,Frequency,Ts);
      MussvInfo.sens = ltipack.frddata(MussvInfo.sens,Frequency,Ts);
      info.MussvInfo = MussvInfo;
      info.MussvBnds = ltipack.frddata(bnds,Frequency,Ts);
   end
   
   % Peak of muLower and muUpper bounds get converted to SM.Upper and Lower
   [pkl,WCIndex] = max(bnds(1,2,:));
   pku = max(bnds(1,1,:));
   stabmarg.UpperBound = inf;
   stabmarg.LowerBound = inf;
   if pkl>0
      stabmarg.UpperBound = 1/pkl;
   end
   if pku>0
      stabmarg.LowerBound = 1/pku;
   end
   stabmarg.DestabilizingFrequency = Frequency(WCIndex); % not in original Units
   
   % Construct worst-case uncertainty
   PertM = mussvunwrap(tinfo.pvec,tinfo.blk);
   if Need.AllBadUncertainValues
      WCU = getWorstCasePerturbation(B,muB,PertM,Frequency);
      DSV = WCU(WCIndex);
      info.BadUncertainValues = WCU;
   else
      DSV = getWorstCasePerturbation(B,muB,PertM(:,:,WCIndex),Frequency(WCIndex));
   end
   
   if Need.Sensitivity
      DisplayFlag = strcmp(Opt.Display,'on');
      gamma = Opt.VaryUncertainty/100;
      idx = find(bnds(1,1,:) >= pku*0.85); %1 x <=npts
      if length(idx) == 1 && npts>1
         if idx == 1
            idx = [idx idx+1];
         else
            idx = [idx-1 idx];
         end
      end
      matg = M11(:,:,idx);
      sfactor = zeros(nblk,1);
      % Options for Sensitivity should be Silent, and take lowerBound down
      % to only 1 try (so multipleTimes = 0).
      SensMuOpt = [Opt.Mussv 's'];
      mloc = find(SensMuOpt=='m');
      if ~isempty(mloc)
         SensMuOpt(mloc+1) = '0';
      end
      for k = 1:nblk
         tmpmat = matg;
         szD = muB.BlkInfo(k).BlockSize;
         affectedRows = muB.BlkInfo(k).ColStart:muB.BlkInfo(k).ColStart+szD(2)-1;
         tmpmat(affectedRows,:,:) = (1+gamma)*tmpmat(affectedRows,:,:);
         name = muB.BlkInfo(k).BlockName;
         if ~strncmpi('rcastreserved',name,13)
            if DisplayFlag
               if k==1
                  t1 = ctrlMsgUtils.message('Robust:analysis:SensitivityDisplay1','ROBUSTSTAB');
                  disp(t1);
               end
               t1 = ctrlMsgUtils.message('Robust:analysis:SensitivityDisplayAll',name);
               disp(['      ' t1 '...']);
            end
            sbnds = mussv(tmpmat,muB.muBlkNby2,SensMuOpt);
            pkus = sbnds(1,1,:); % 1-by-1-by-(1,2 or more)
            oldpeak = bnds(1,1,idx);
            oldpeak = oldpeak(:); %col vector
            pkus = pkus(:);
            if any(oldpeak==0 & pkus>0)
               sfactor(k) = inf;
            elseif all(oldpeak==0 & pkus==0)
               sfactor(k) = 1;
            else
               keepgz = find(oldpeak>0);
               sfactor(k) = max(abs( (pkus(keepgz)-oldpeak(keepgz)) ./ (gamma*oldpeak(keepgz)) ));
            end
         end
      end
      SensValues = round(100*sfactor);
      info.Sensitivity = cell2struct(num2cell(SensValues),BlockName,1);
      SensText = rctutil.makeSensMessage(info.Sensitivity,1); % CHAR array for Report
   else
      info.Sensitivity = cell2struct(num2cell(NaN([nblk 1])),BlockName,1);
      SensText = '';
   end
end