function plotgauss(C)

% make a contour plot of a bivariate Gaussian with covariance matrix C

% Chris Williams, November 2009


[nr, nc] = size(C);

% check size of C
if (nr ~= 2 || nc ~=2)
  error('input matrix must be 2x2')
end

% set up grid 

s1 = sqrt(C(1,1));  % std deviation on axis 1
s2 = sqrt(C(2,2));  % std deviation on axis 2

x = linspace(-2*s1,2*s1);	% location of points at which x is calclated
y = linspace(-2*s2,2*s2);	% location of points at which y is calclated

[X, Y] = meshgrid(x,y); % matrices used for plotting


% compute quadratic form

P = inv(C);		   % the inverse covariance matrix

Z = P(1,1)*X.^2 + 2*P(1,2)*X.*Y + P(2,2)*Y.^2;  % calculation of the quadratic
						% form

% compute Gaussian and make contour plot

ncontour = 20;  % number of contours

G = exp(- 0.5*Z);		% Note: Gaussian is not normalized, but this
				% doesn't matter for plotting purposes
figure(1)
contour(x,y,G,ncontour)	% contour plot
axis('square')			% command to set axes square, so circles
				% look like circles