function R=are_schur(A,B,Q,P)
% ARE_SCHUR Solve the continuous Algebraic Riccati Equation (ARE)
% with the block real schur decomposition of the Hamiltonian matrix.
% SYNTAX
% R=are_schur(A,B,Q,P)
% ON INPUT
% A,B - system matrices
% Q   - state weighting matrix
% P   - control weighting matrix
% ON OUTPUT
% R   - The positive definite solution to the ARE, A'R+RA-RBinv(P)B'R+Q=0.
% NON STANDARD MATLAB FUNCTIONS CALLED:
% BLKRSCH2, CSCHUR2
% Written by David Di Ruscio

F=[ A -B*pinv(P)*B';-Q -A'];           % Hamilton matrix
[U,T]=blkrsch2(F,1);                   % stable eigenvalues in upper
                                       % n x n submatrix T_11.
n=length(A);
U11=U(1:n,1:n);
U21=U(n+1:2*n,1:n);
R=U21*pinv(U11);

