function [Value,RespObj] = calcopt(Property,RespObj,varargin)
%CALCOPT calculates Plot Options for Impulse 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 = CALCOPT(Property,RespObj,Index) calculates the values for
%   only the systems in RespObj for the indices in Index.
%
%   [Value,RespObj] = CALCOPT(Property,RespObj) also returns the 
%   Response Object with the Property value set to the calculated
%   value.

%   Copyright 1986-2000 The MathWorks, Inc. 
%   Karen Gondoly 1-28-97.
%   $Revision: 1.10 $  $Date: 2000/06/14 20:10:51 $

RespProps = get(RespObj,{'ResponseHandles','SystemNames'});
[ResponseHandles,SysNames] = deal(RespProps{:});

if nargin>2,
   Ind = varargin{1};
else
   Ind = [1:length(ResponseHandles)];   
end

switch lower(Property),
case 'peakresponsevalue',
   Value = struct('System',cell(length(Ind),1),'Time',[],'Peak',[]);
case 'settlingtimevalue',
   Value = struct('System',cell(length(Ind),1),'SettlingTime',[],'Amplitude',[]);
   SetTimeLims = RespObj.SettlingTimeThreshold;
end

% Assign uptodate system names to Value.System field
% (In the LTI Viewer, the names are not available when RespObj is created)
[Value(Ind).System] = deal(SysNames{Ind});

% Compute required characteristic value
for ctInd = 1:length(Ind);
   AllHandles = ResponseHandles{Ind(ctInd)}; % Handles for a particular system
   Xvals = zeros([size(AllHandles),size(AllHandles{1})]);
   Yvals = zeros(size(Xvals));
   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');
            Y = get(ResponseLine,'Ydata');
            absY=abs(Y);
            
            switch lower(Property),
            case 'peakresponsevalue',
               [garb,indMax]=max(absY);
               Xvals(ctrow,ctcol,ctModel) = X(indMax);
               Yvals(ctrow,ctcol,ctModel) = Y(indMax);
               
            case 'settlingtimevalue',
               % RE: will give wrong answer for unstable models!
               MaxError = max(abs(Y));
               NotSettled = find(abs(Y)>SetTimeLims*MaxError);
               if isempty(NotSettled),
                  % Direct feedthrough
                  Ts = 0;
                  As = Y(1);
               elseif NotSettled(end)==length(Y)
                  % Not Settled yet
                  Ts = NaN;
                  As = 0;
               else   
                  idx = NotSettled(end);
                  if Y(idx)<0,
                     As = -SetTimeLims*MaxError;
                  else
                     As = SetTimeLims*MaxError;
                  end
                  Ts = X(idx) + (X(idx+1)-X(idx))/(Y(idx+1)-Y(idx)) * (As-Y(idx));
               end 
               Xvals(ctrow,ctcol,ctModel) = Ts;
               Yvals(ctrow,ctcol,ctModel) = As;
            end % switch Property
         end % for ctModel
      end % for ctcol
   end % for ctrow
   
   %---Store the calculated values into the structured array
   switch lower(Property),
   case 'peakresponsevalue',
      Value(ctInd).Time = Xvals;
      Value(ctInd).Peak = Yvals;
      RespObj.PeakResponseValue = Value;
   case 'settlingtimevalue',
      Value(ctInd).SettlingTime = Xvals;
      Value(ctInd).Amplitude = Yvals;
      RespObj.SettlingTimeValue = Value;
   end
   
end % for ctInd

