function [Value,RespObj] = calcopt(Property,RespObj)
%CALCOPT calculates Plot Options for Nichols Response Plots
%
%   Value = CALCOPT(Property,RespObj) calculates the values of the
%   Plot Options specified by Property for all systems in the Response
%   Object, RespObj. The values are returned as the cell array, Value,
%   where each cell represents the results for each system in RespObj.
%
%   [Value,RespObj] = CALCOPT(Property,RespObj) also returns the 
%   Response Object with the Property value set to the calculated
%   value.

%   Authors: Karen Gondoly, 1-28-97
%   Revised: Adam DiVergilio, 11-99
%   Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.18 $  $Date: 2000/06/14 20:09:38 $

ResponseHandles = get(RespObj,'ResponseHandles');
SysNames = get(RespObj,'SystemNames');
Nsys = length(SysNames);

%---Turn warnings off, temporarily, to avoid MARGIN warnings
WarnStr = warning; warning off;

switch lower(Property)
case 'peakresponsevalue'
   Value = struct('System',cell(Nsys,1),'Frequency',[],'Peak',[],'PhaseAtPeak',[]);
   for ctSys = 1:Nsys
      Value(ctSys).System = SysNames{ctSys};
      AllHandles = ResponseHandles{ctSys};  % Handles for a particular system
      Freqs = ones([size(AllHandles),size(AllHandles{1})]);
      Peaks = ones(size(Freqs));
      PeakPhase = ones(size(Freqs));
      for ctrow = 1:size(AllHandles,1)
         for ctcol = 1:size(AllHandles,2)
            Handles = AllHandles{ctrow,ctcol};
            %---Go through each model
            for ctModel = 1:prod(size(Handles))
               ResponseLine = findobj(Handles{ctModel},'Tag','LTIresponseLines');
               X = get(ResponseLine,'Xdata');  % phase
               Y = get(ResponseLine,'Ydata');  % mag
               [Peaks(ctrow,ctcol,ctModel),indMax] = max(Y);
               Freqs(ctrow,ctcol,ctModel) = RespObj.Frequency{ctSys}{ctModel}(indMax);
               PeakPhase(ctrow,ctcol,ctModel) = X(indMax);
            end
         end
      end
      Value(ctSys).Frequency = Freqs;
      Value(ctSys).Peak = Peaks;
      Value(ctSys).PhaseAtPeak = PeakPhase;
   end
   RespObj.PeakResponseValue = Value;

case 'stabilitymarginvalue'
   Value = struct('System',cell(Nsys,1),...
      'Stable',[],'GainMargin',[],'GMFrequency',[],'PhaseMargin',[],'PMFrequency',[]);
   LTIviewerFig = get(RespObj,'Parent');
   if strcmp(get(LTIviewerFig,'Tag'),'ResponseGUI')
      ViewerObj = get(LTIviewerFig,'UserData');
      Systems = get(ViewerObj,'Systems');
   else
      Nsys = 0;
   end
   for ctSys = 1:Nsys
      Value(ctSys).System = SysNames{ctSys};
      if ~isa(Systems{ctSys},'lti')
         Stable=[]; Gm=[]; Pm=[]; Wcg=[]; Wcp=[];
      else
         [Gm,Pm,Wcg,Wcp,Stable] = margin(Systems{ctSys});
      end
      Value(ctSys).Stable = Stable;
      Value(ctSys).GainMargin  = unitconv(Gm,'abs',     RespObj.MagnitudeUnits);
      Value(ctSys).PhaseMargin = unitconv(Pm,'deg',     RespObj.PhaseUnits);
      Value(ctSys).GMFrequency = unitconv(Wcg,'rad/sec',RespObj.FrequencyUnits);
      Value(ctSys).PMFrequency = unitconv(Wcp,'rad/sec',RespObj.FrequencyUnits);
   end
   RespObj.StabilityMarginValue = Value;

case 'allstabilitymarginvalue'
   Value = struct('System',cell(Nsys,1),...
      'Ts',cell(Nsys,1),'Margins',cell(Nsys,1));
   LTIviewerFig = get(RespObj,'Parent');
   if strcmp(get(LTIviewerFig,'Tag'),'ResponseGUI')
      ViewerObj = get(LTIviewerFig,'UserData');
      Systems = get(ViewerObj,'Systems');
   else
      Nsys = 0;
   end
   for ctSys = 1:Nsys
      Value(ctSys).System = SysNames{ctSys};
      if ~isa(Systems{ctSys},'lti')
         Ts=[]; S=[];
      else
         Ts = Systems{ctSys}.Ts;
         S = allmargin(Systems{ctSys});
         syssize = size(Systems{ctSys});
         nmodels = prod(syssize(3:end));
         for n=1:nmodels
            S(n).GainMargin  = unitconv(S(n).GainMargin, 'abs',    RespObj.MagnitudeUnits);
            S(n).PhaseMargin = unitconv(S(n).PhaseMargin,'deg',    RespObj.PhaseUnits);
            S(n).GMFrequency = unitconv(S(n).GMFrequency,'rad/sec',RespObj.FrequencyUnits);
            S(n).PMFrequency = unitconv(S(n).PMFrequency,'rad/sec',RespObj.FrequencyUnits);
            S(n).DMFrequency = unitconv(S(n).DMFrequency,'rad/sec',RespObj.FrequencyUnits);
            if length(S(n).DelayMargin)>length(S(n).PhaseMargin)
               if isinf(S(n).DMFrequency(end))
                  S(n).DMResponse = evalfr(Systems{ctSys}(:,:,n),inf);
               else
                  S(n).DMResponse = evalfr(Systems{ctSys}(:,:,n),-1);
               end
            else
               S(n).DMResponse = [];
            end
         end
      end
      Value(ctSys).Ts = Ts;
      Value(ctSys).Margins = S;
   end
   RespObj.AllStabilityMarginValue = Value;

end

%---Restore warning state
warning(WarnStr);
