function [mu, Sigma] = add_ev_to_gaussian(muX, SigmaX, X, Y, dom, y, sizes)
% ADD_EV_TO_GAUSSIAN Compute Pr(X,Y=y) where X and Y are jointly Gaussian.
% [mu, Sigma] = add_ev_to_gaussian(muX, SigmaX, X, Y, domain, y, sizes)
% 
% sizes(i) is the size of the i'th block in domain
% mu(X) = muX, mu(Y) = y
% Sigma(X,X) = SigmaX, Sigma(X,Y) = Sigma(Y,X) = Sigma(Y,Y) = 0

n = sum(sizes);
mu = zeros(n, 1);
Sigma = zeros(n,n);

if ~isempty(X)
  i = find_equiv_posns(X, dom);
  b = block(i, sizes);
  mu(b) = muX;
  Sigma(b,b) = SigmaX;
end

if ~isempty(Y)
  i = find_equiv_posns(Y, dom);
  b = block(i, sizes);
  mu(b) = y(:);
end


