function [Ts,Tf] = timegrid(a,b,c,x0,MaxDelay,Nsmax,Tf)
%TIMEGRID  Estimates step size and time horizon for continuous-time responses.
%
%   [TS,TF] = TIMEGRID(A,B,C,X0,MAXDELAY,NSMAX) computes an adequate fixed   
%   step size TS for simulating the time response of the affine ODE
%       dx/dt = Ax + B       x(0) = X0
%          y  = Cx
%   TIMEGRID also returns an estimate TF of how long the response takes
%   to settle. For unstable systems, TF is the suggested simulation horizon.  
%   B and X0 can be matrices, in which case TIMEGRID selects adequate TS 
%   and TF for the set of simulations with step input B(:,1),B(:,2),... 
%   and initial conditions X0(:,1), X0(:,2),... NSMAX is the maximum number
%   of time steps allowed during simulation.
%
%   TS = TIMEGRID(A,B,C,X0,MAXDELAY,TF,NSMAX) estimates the step size TS 
%   given the simulation horizon TF in seconds.
%
%   TIMEGRID attempts to produce a simulation time TF that ensures the
%   output responses have decayed to approximately 1% of their initial
%   or peak values. 

%   Author: P. Gahinet
%   Copyright 1986-2012 The MathWorks, Inc.

% Parameters
ZeroTol = 1e-8;
Nsmin = 100;  % min number of time steps

% Note: Currently coded assuming b=[] or x0=[]
nx = size(a,1);
bx0 = [b,x0];
Tfmin = 2*MaxDelay;
StepFlag = ~isempty(b);
NoTf = (nargin<7 || isempty(Tf));

% Enforce consistent scaling of output channels
if StepFlag && size(c,1)>1 && rcond(a)>ZeroTol
   cs = max(abs(c*(a\bx0)),[],2);
   cs(cs==0) = 1;
   c = lrscale(c,1./cs,[]);
end

% Compute poles and detect perturbed integrators
[p,v,w,pcond] = localSpectralDecomp(a);
wp = abs(p);
rp = real(p);
rpmax = max(rp);  % most unstable mode
isSimpleMode = (pcond<1e8);

% Compute DC contribution of each mode
% Note: Contribution of multiple eigenvalues can't be estimated (Contrib=0)
beta = w'*bx0;
gamma = c*v;
Contrib = zeros(nx,1);
for ct=1:nx
   if isSimpleMode(ct)
      Contrib(ct) = norm(gamma(:,ct)) * norm(beta(ct,:)) * pcond(ct);
   end
end
Contrib(imag(p)~=0,:) = 2*Contrib(imag(p)~=0,:);  % complex pairs

% Adjust contribution for step response
if StepFlag
   % Note: Near-cancelling modes near s=0 have 0/0 contribution. Force
   % contribution to zero when (c*vj)*(wj'*b) is small (see cancel01)
   Contrib(Contrib<ZeroTol*max(Contrib) & wp<ZeroTol*min(1,max(wp))) = 0;
   Contrib = Contrib ./ (eps+wp);  % prevent Inf's
end

% Flag integral action
IntegralAction = any(wp==0 & Contrib>0.1*max(Contrib));
Contrib(wp==0) = 0;  % ignore integrators when picking Tf,Ts

% Quick exit if pure gain or integrator
if all(Contrib==0)
   if NoTf  
      Tf = max(1,Tfmin);  
   end
   Ts = Tf/Nsmin;
   return
end

% Dominant modes are non-integrating modes contributing at least 10%
iDom = find(wp>0 & (Contrib>0.1*max(Contrib) | ~isSimpleMode));
wp = wp(iDom,:);
rp = rp(iDom,:);

% Compute adequate Ts,Tf for dominant modes
% Tf selection: Treat response as exp(j*wp) or exp(rp*t)
isJW = (abs(rp)<ZeroTol*wp);                          % pure imaginary
Tfp(isJW,:) = (10*pi)./wp(isJW,:);                    % ~5 periods
Tfp(~isJW,:) = log(100)./(1e3*eps+abs(rp(~isJW,:)));  % factor 100 change in exp(r*t)
% Ts selection: Sample based on time constants and oscillation periods
Tsp = Tfp/50;
isRes = (wp.*Tfp>30);  % more than 5 oscillations periods in [0,Tf]
Tsp(isRes,:) = min(Tsp(isRes,:),(pi/10)./wp(isRes,:));   % ~20 points/period

% Pick final time if not supplied
if NoTf
   Tf = max(Tfp);
   % Simulate long enough to show integral action and just enough for
   % most unstable mode to diverge
   if IntegralAction
      Tf = 10*Tf;
   end
   if rpmax>ZeroTol
      Tf = min(Tf,log(1e8)/rpmax);
   end
   % Enforce minimum duration
   Tf = max(Tfmin,Tf);
end

% Pick sample time
Ts = min(Tsp);
% Enforce minimum number of steps
Ts = min(Ts,Tf/Nsmin);

% Out of memory safeguards
if NoTf
   % Unspecified final time: protect against ultrafast dynamics or very HF modes
   Ts = max(Ts,Tf/Nsmax);
else
   % Specified final time: protect against final time much larger than
   % system's time constant (can happen in FILLRESP, see g258770)
   Ns = min(round(Tf/Ts),Nsmax);
   Ts = Tf/Ns;  % to ensure simulation actually stops at Tf
end

%[Ts Tf]

%-------------------------------

function [p,v,w,pcond] = localSpectralDecomp(a)
% Compute eigenstructure and detect pseudo-integrators
nx = size(a,1);
if nx==0
   p = zeros(0,1);  v = [];  w = [];   pcond = zeros(0,1);
else
   % Compute left and right eigenvectors
   [s,perm,a] = balance(a);
   [v,p,w] = eig(a,'vector');
   pcond = zeros(nx,1);
   for ct=1:nx
      pcond(ct) = 1 / abs(w(:,ct)' * v(:,ct));
   end
   % Zero out eigenvalues within rounding errors of s=0
   % g1342285,g1451123: BALANCING needed to avoid zeroing most eigs 
   %                    when ||A||>>rho(A)
   ztol = min(1e10,max(1e4,pcond)) * eps * norm(a,1);
   p(abs(p)<ztol) = 0;
   % Undo balancing
   w = lrscale(w(perm,:),1./s,[]);
   v = lrscale(v(perm,:),s,[]);
end
