function  [sys, x0, str, ts]  = rlsests(t,x,u,flag,nstates,lambda,dt)
%RLSESTS S-function to perform system identification.
%
%   This M-file is designed to be used in a Simulink S-function block.
%   It performs parameter estimation using the Recursive Least Squares
%   Parameter Estimation Algorithm with Exponential Data Weighting
%   
%   The input arguments are
%
%     nstates:        the number of states in the states vector
%     lambda:         the exponential data weighting factor
%     dt:             how often to sample points (secs)
%
%   The RLS estimator is defined by the following equations:
%
%   	                        1      P(k-2) * phi(k-1) * [y(k) - phi(k-1)'theta(k-1)]
%     theta[k] = theta[k-1] + ------ * -------------------------------------------------
%   	                      lambda         lambda + phi(k-1)' * P(k-2) * phi(k-1)
%
%   	         1       P(k-2) * phi(k-1) * phi(k-1)' * P(k-2)
%     P(k-1) = ------ * ----------------------------------------
%   	       lambda    lambda + phi(k-1)' * P(k-2) * phi(k-1)
%
%   where:
%
%     theta:		the parameter estimates
%     phi:		the state vector
%     P:		the covariance matrix
%     lambda:		the exponential data weighting factor
%
%   See also SFUNTMPL., "Adaptive Filtering, Prediction, and Control",
%   G. C. Goodwin & K. S. Sin.

%   Copyright (c) 1990-1998 by The MathWorks, Inc. All Rights Reserved.
%   $Revision: 1.15 $
%   Rick Spada 6-17-92.


if abs(flag) == 2	% flag = 2 --> real time hit
  
  % sample hit, return the next discrete states, which are the
  % next parameter estimates

  theta = x(1:nstates);					% parameter estimates
  P = zeros(nstates,nstates);				% get covariance matrix
  P(:) = x(nstates+1:nstates+nstates*nstates);
  yk = u(nstates + 1);					% system output
  phi = u(1:nstates);					% state vector
  est_err = yk - phi' * theta;				% estimation error
  den = lambda + phi' * P * phi;			% lambda + phi' * P * phi
  theta_new = theta + P * phi * (est_err / den);	% new parameter estimates
  Pnew = (P - P * phi * phi' * P / den) / lambda;	% new covariance
  sys = [theta_new', Pnew(:)']';			% return them
  
elseif flag == 4	% flag = 4 --> Return next sample hit

  sys = [];

elseif flag == 0	% flag = 0 --> Return sizes of parameters and initial conditions
  
  sys(1) = 0;				% 0 continuous states
  sys(2) = nstates+nstates*nstates;	% enough discrete states to hold the estimates
  					% and the covariance matrix
  sys(3) = nstates;			% nstate estimate outputs
  sys(4) = nstates+1;			% nstate+1 (regression vector + system output)
  					% inputs
  sys(5) = 0;				% 0 roots
  sys(6) = 0;				% no direct feedthrough
  sys(7) = 1;				% 1 sample time

  % initialize the covariance matrix and initial estimates
  P = eye(nstates, nstates) * 1e6;
  x0 = [ones(nstates,1)', P(:)']';
  ts = [dt, 0];
  
elseif flag == 3	% flag = 3 --> Return outputs, only at sample hits
  
    sys = x(1:nstates);
    sys = sys(:); 
else

  sys = [];
  
end
