% Ricatti evaluation file
% =======================

function result = RicattiEval(x0,p,Y)

global hetds

Q0U = hetds.Q0;
Q1S = hetds.Q1;
A = Het_mjac(x0,p);
result = [];

if ~isempty(YU)
    % Riccati blocks from unstable eigenspace
    [U11, U12, UE21, U22] = RicattiCoeff(Q0U,A,hetds.nu);
    tmp = (U22*YU - YU*U11 + UE21 - YU*U12*YU)';
    for i=1:hetds.nu
        result(i,1) = tmp(:,i);
        %result(end+1:end+hetds.npos,1) = tmp(:,i);
    end
end

if ~isempty(YS)
    % Riccati blocks from stable eigenspace
    [S11, S12, SE21, S22] = RicattiCoeff(Q0S,A,hetds.ns);
    tmp = (S22*YS - YS*S11 + SE21 - YS*S12*YS)';
    for i=1:hetds.ns
        result(i,1) = tmp(:,i);
    end
end