function [Value,RespObj] = calcopt(Property,RespObj,varargin)
%CALCOPT calculates Plot Options for Step 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.

% $Revision: 1.13 $
%   Copyright 1986-2000 The MathWorks, Inc. 
%   Karen Gondoly 1-28-97.


SteadyStateVal = RespObj.SteadyStateValue;
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 'risetimevalue',
   Value = struct('System',cell(length(Ind),1),...
      'StartTime',cell(length(Ind),1),...
      'RiseTime',cell(length(Ind),1),...
      'Amplitude',cell(length(Ind),1));
   RiseTimeLims = RespObj.RiseTimeLimits;
      
case 'settlingtimevalue',
   Value = struct('System',cell(length(Ind),1),'SettlingTime',cell(length(Ind),1),...
      'Amplitude',cell(length(Ind),1));
   SetTimeLims = RespObj.SettlingTimeThreshold;
   
case 'steadystatevalue',
   Value = SteadyStateVal;
   
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});

% Exit if steady state value required
% RE: CALCOPT only sets the system names for 'steadystatevalue'
if strcmpi(Property,'steadystatevalue')
   return
end

% Compute required characteristic values
for ctInd = 1:length(Ind);
   AllHandles = ResponseHandles{Ind(ctInd)};    % Handles for a particular system
   K = SteadyStateVal(Ind(ctInd)).Amplitude;    % DC gain (may be infinite)
   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)),
            DC = K(ctrow,ctcol,ctModel);   % DC gain for current model and I/O pair
            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 'risetimevalue',
               if isinf(DC),
                  % Infinite DC gain
                  Xvals(ctrow,ctcol,ctModel) = Inf;
                  Yvals(ctrow,ctcol,ctModel) = Y(end);
                  Startvals(ctrow,ctcol,ctModel) = Inf;
               else
                  % Finite DC
                  Y0 = Y(1);  % initial value
                  Y1 = Y0 + RiseTimeLims(1)*(DC-Y0);  % lower threshold line
                  Y2 = Y0 + RiseTimeLims(2)*(DC-Y0);  % upper threshold line
                  % Find first crossing for each threshold line
                  Npts = length(Y);
                  IndStart = find(sign(Y(1:Npts-1)-Y1).*sign(Y(2:Npts)-Y1)<=0);
                  IndEnd = find(sign(Y(1:Npts-1)-Y2).*sign(Y(2:Npts)-Y2)<=0);
                  % Get time of first crossing
                  if isempty(IndStart),
                     % Has not yet reached RiseTimeLims(1) level
                     X1 = NaN;
                  else
                     idx = IndStart(1);
                     if Y(idx)==Y(idx+1),
                        X1 = X(idx);
                     else
                        X1 = X(idx) + (X(idx+1)-X(idx))/(Y(idx+1)-Y(idx)) * (Y1-Y(idx));
                     end
                  end
                  % Get time of second crossing
                  if isempty(IndEnd),
                     % Has not yet reached RiseTimeLims(2) level
                     Tr = NaN;
                  else
                     idx = IndEnd(1);
                     if Y(idx)==Y(idx+1),
                        X2 = X(idx);
                     else
                        X2 = X(idx) + (X(idx+1)-X(idx))/(Y(idx+1)-Y(idx)) * (Y2-Y(idx));
                    end
                    Tr = X2-X1;
                end
                Xvals(ctrow,ctcol,ctModel) = Tr;
                Yvals(ctrow,ctcol,ctModel) = Y2;
                Startvals(ctrow,ctcol,ctModel) = X1;
            end
               
            case 'settlingtimevalue',
               if isinf(DC),
                  % No finite settling time
                  Ts = Inf;
                  As = Y(end);
               else
                  MaxError = max(abs(Y-DC));
                  NotSettled = find(abs(Y-DC)>SetTimeLims*MaxError);
                  if isempty(NotSettled),
                     % Direct feedthrough
                     Ts = 0;
                     As = Y(1);
                  elseif NotSettled(end)==length(Y)
                     % Not Settled yet
                     Ts = NaN;
                     As = DC;
                  else   
                     idx = NotSettled(end);
                     if Y(idx)<DC,
                        As = DC-SetTimeLims*MaxError;
                     else
                        As = DC+SetTimeLims*MaxError;
                     end
                     Ts = X(idx) + (X(idx+1)-X(idx))/(Y(idx+1)-Y(idx)) * (As-Y(idx));
                  end 
               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 'risetimevalue',
      Value(ctInd).StartTime = Startvals;
      Value(ctInd).RiseTime = Xvals;
      Value(ctInd).Amplitude = Yvals;
      RespObj.RiseTimeValue = Value;
   case 'settlingtimevalue',
      Value(ctInd).SettlingTime = Xvals;
      Value(ctInd).Amplitude = Yvals;
      RespObj.SettlingTimeValue = Value;
   end
   
end % for ctInd

% end ../@stepplot/calcopt.m
