function [lyapunov,t] = analyticallyapunov(ode,tspan,yo,odeopts,nsigma,tstep,varargin)
% ANALYTICALLYAPUNOV Estimate Lyapunov exponents given dynamic equations.
%   LYAPUNOV = ANALYTICALLYAPUNOV(ODE,TSPAN,YO) with TSPAN = [T0 TFINAL] 
%   calculates the Lyapunov exponents for a continuous system by integrating the 
%   variational equations (jacobian) along the solution from inital condition, 
%   YO. Default renormalization occurs 100 times over the interval TSPAN. The 
%   continuous system's and variational equation's dervitives, Y'=f(t,Y), are 
%   defined in the function specified by the function handle ODE. The solver is 
%   ODE45.
% 
%   For a system of dimension n, Y'=f(t,Y) should include a total of (n+n^2) 
%   equations. n equations for the evolution of the trajectory from Y0, and n^2 
%   equations describing the evolution of the variational equations. (see
%   example of Lorenz equations, ./EXAMPLES/DLORENZ.M)
%
%   LYAPUNOV = ANALYTICALLYAPUNOV(ODE,TSPAN,YO,ODEOPTS) pass options created 
%   with ODESET to the ode45 solver.
%
%   LYAPUNOV = ANALYTICALLYAPUNOV(ODE,TSPAN,YO,ODEOPTS,NSIGMA) where NSIGMA is a 
%   n x 1 column vector that is normal to the Poincare section that halted
%   integration. This is useful for hybrid systems that have a discrete event 
%   that is applied at intersection with the Poincare section. Pass an empty 
%   argument, [], to not use this option.
%
%   LYAPUNOV = ANALYTICALLYAPUNOV(ODE,TSPAN,YO,ODEOPTS,NSIGMA,TSTEP) where TSTEP
%   is the time step to perform renormalization. Pass an empty argument,[], to 
%   use the default of 100 renormalizations over TSPAN.
%
%   LYAPUNOV = ANALYTICALLYAPUNOV(ODE,TSPAN,YO,ODEOPTS,NSIGMA,TSTEP,...) any 
%   additional arguments will be passed directly to the ODE45 solver.
%
%   [LYAPUNOV,T] = ANALYTICALLYAPUNOV(...) Outputs corresponding time of
%   evolution.
%
%   This function calls ANALYTICALFLOQUET to calculate the state-transition 
%   matrix. ANALYTICALFLOQUET provides a warning if the initial condition does 
%   not agree with the final state after integration, which is normally the case 
%   because Lyapunov exponents are only intersting for aperiodic data. By 
%   default, this warning is disabled when this function is called and 
%   re-enabled after exection. See: warning('off', 'analyticalfloquet:refsol')         
% 
%   Input
%       ODE      | 1 x 1 | function handle | e.g. @dynamics
%       TSPAN    | 1 x 2 | double          | simulation time [TO TFINAL]
%       YO       | n x 1 | double          | initial state
%       ODEOPTS  | 1 x 1 | struct          | ode options, use: odeset
%       TSTEP    | 1 x 1 | double          | time for renormalization
%       NSIGMA   | n x 1 | double          | vector normal to Poincare section
%       VARARGIN |   ?   |                 | passed directly to odesolver       
%
%   Output
%       LYAPUNOV | m x n | double | Estimates of Lyapunov exponents
%       T        | m x 1 | double | Time corresponding to Lyapunov
%   
%   Example        
%       lorenzex; % Run the example script
%
%   Reference
%     [1] Wolf, A, Swift, JB, Swinney, HL, Vastano, JA, "Determining Lyapunov 
%         exponents from a time series," Physica D: Nonlinear Phenomena, Vol. 
%         16, No. 3, pp. 285-317, 1985. 
%
%     [2] Hiskens, IA, "Stability of hybrid system limit cycles: application to 
%         the compass gait biped robot," In 40th IEEE Conference on Decision and 
%         Control, 2001, pp. 774-779.
%
%   See also LORENZEX, DLORENZ, ANALYTICALFLOQUET.

% Created by: Jesse Norris
%             Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%             http://jesse.musculoskeletaldynamics.com
% Last Revision: 3 October 2007

% Turn off warning for non-period solutions
  warning('off', 'analyticalfloquet:refsol')
  
% Handle the optional input arguments dSquared and nNearest
  if ~exist('odeopts','var') || isempty(odeopts); odeopts = odeset; end;
  if ~exist('tstep','var')   || isempty(tstep);   tstep   = diff(tspan)/100; end;
  if ~exist('nsigma','var'); nsigma = []; end;
  
% Initialize parameters
   n    = length(yo);      % number of states
    
   Q = eye(n);              % Initial rotation matrix
   cumulative = zeros(1,n); % Cumulative tracking of disturbances after QR normalizaiton   
   t          = (tspan(1):tstep:tspan(end))'; % Time
   lyapunov   = zeros(length(t)-1,n);         % Lyapunov exponents at each step
   
  
% Main loop    
   for i = 1:length(t)-1       
       % Estimate state-transition matrix phi
         [junk,phi,yo] = analyticalfloquet(ode,t([i i+1]),yo(:),odeopts,nsigma,varargin{:});                       
         
       % Actual estimation of Lyapunov exponents  
       % Note that phi is multiplied by the previous Q since the function
       % analyticalfloquet always begins from the identity matrix.
         [Q,R] = qr(phi*Q); % Gram-Schmidt orthonormalization Q=orthonorm vector, R=magnitude
         cumulative = cumulative + log(abs(diag(R)))';  % vector length growth = abs(diag(R))       
         lyapunov(i,:) = cumulative/diff(t([1 i+1]));   % Lyapunov exponents      
   end;  
   
% Truncate time since you don't have an estimate of Lyapunov spectrum at t(end)
  t = t(1:end-1);
   
% Turn on warning for non-period solutions 
  warning('on', 'analyticalfloquet:refsol')

end