function[a,b,g]=mult(a1,b1,g1,a2,b2,g2)
%function[a,b,g]=mult(a1,b1,g1,a2,b2,g2)
%  Multiplies two symplectic matrices.  It creates E2*E1 where Ei is
%  represented in the form
%
%            |                                    |
%            |   ai^(-1)         ai^(-1)*bi       |
%        Ei =|                                    |         (&)
%            |  gi*ai^(-1)    ai'+gi*ai^(-1)*bi   |
%            |                                    |
%
%  It creates E=E2*E1 where E is represented in the same form as Ei.
%  The algorithm is
%
%        a = a2*[(I + b1*g2)^(-1)]*g1
%        g = g1 + a1'*g2*[(I+b1*g2)^(-1)]*a1
%        b = b2 + a2*[(I + b1*g2)^(-1)]*b1*a2'
%
%  where a, g, and b are the objects defining E in the counterpart to (&).
%  See Anderson and Moore, Optimal Filtering, page 160 for the definition
%  of symplectic matrices and for a demonstration that they can be represented
%  in this way.
%  The matrices (a1,a2), (b1,b2), and (g1,g2), respectively, must each be
%  of the same dimensions.
ss=max(size(a1));
v=eye(ss);
a=a1*((v+b2*g1)\a2);
g=g2+a2'*g1*((v+b2*g1)\a2);
b=b1+a1*((v+b2*g1)\b2*a1');
