function [K,S]=kfricc(A,C,Q,R,W);
%KFRICC  Iteration on Riccati to solve for Kalman filter.
%        [K,S]=kfricc(A,C,Q,R,W) calculates the Kalman gain, K, and
%        the stationary covariance matrix, S, obtained from the 
%        Kalman filter with the following system:
%
%	  x[t+1] = Ax[t] + Bu[t] + w1[t+1]
%           y[t] = Cx[t] + Du[t] + w2[t]
%                                                          (1)
%         E [w1(t+1)] [w1(t+1)]' =  [Q  W;
%           [ w2(t) ] [ w2(t) ]      W' R]
%
%        where x is the mx1 vector of states, u is the nx1 vector of 
%        controls, y is the px1 vector of observables, A is mxm, B is
%        mxn, C is pxm, Q is mxm, R is pxp, and W is mxp.  The 
%        innovations representation of system (1) is then
%         _               _                     _
%         x[t+1] = (A-K*C)x[t] + Bu[t] +Ky[t],  x[t]=E[ x[t] | y[t],y[t-1]..]
%                         _   
%           y[t] = Cx[t] +Du[t] +a[t]                            
%                                                       _            _
%        where K = (A*S*C'+W)/(R+C*S*C'),  S= E[(x(t)-x(t))*(x(t)-x(t))']
%
m=max(size(A));
[rc,cc]=size(C);
if nargin<5; W=zeros(m,rc); end;
S=.01*eye(m);
k0=(A*S*C'+W)/(R+C*S*C');
dd=1;
it=1;
maxit=5000;
tol=1e-8;       % reset tolerance if necessary 
while (dd> tol & it<= maxit);
  S =A*S*A' + Q -(A*S*C'+W)*k0';
  k1=(A*S*C'+W)/(R+C*S*C');
  dd=max(max(abs( (k1-k0)./(1+abs(k0)) )));
  it=it+1;
  k0=k1;
end;
K=k1;
if it>=maxit; disp('WARNING: Iteration limit of 5000 reached in DRICC'); end;
