function plotquad(v1, v2, rho)

% make a contour plot  of the quadratic form that defines a bivariate
% Gaussian, and a surface plot of the corresponding Gaussian

% Chris Williams, October 1999

% plot quadratic using the meshgrid commands

s1 = sqrt(v1);  % std deviation on axis 1
s2 = sqrt(v2);  % std deviation on axis 2

if (abs(rho) >= 1.0)
  disp('error: rho must lie between -1 and 1');
  return
end
cov12 = rho*s1*s2;  % calculation of the covariance 

C = [v1 cov12; cov12 v2];  % the covariance matrix
A = inv(C);		   % the inverse covariance matrix

x = -2*s1:0.1:2*s1;	% location of points at which x is calclated
y = -2*s2:0.1:2*s2;	% location of points at which y is calclated

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

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

figure(1)
contour(x,y,Z, [1 1], 'b-')	% contour plot
axis('square')			% command to set axes square, so circles
				% look like circles

G = exp(- 0.5*Z);		% Note: Gaussian is not normalized, but this
				% doesn't matter for plotting purposes
figure(2)
surf(x,y,G)