#Code adapted from Craig John's code bivnorm.R
## Enter the Sigma matrix 
SG <- matrix(c(4,-3,-3,9),2,2);

## SGI is the inverse of Sigma (variance-covariance matrix)
SGI <- solve(SG);

## The mean vector
MU = c(1,0)

## The normalizing constant for the bivariate normal denisty
kk = 2*pi*sqrt(SG[1]*SG[4]-SG[3]^2)

## A set of x and y locations at which to evaluate the density
xx = seq(-1,1,,100)*sqrt(SG[1,1])*2.5 + MU[1];
yy = seq(-1,1,,80)*sqrt(SG[2,2])*2.5 + MU[2];


zz <- matrix(0,100,80);
for( i in 1:100){
for( j in 1:80){
  dum     <- c(xx[i]-MU[1],yy[j]-MU[2]);       # (location - mean)
  zz[i,j] <- exp(-.5*(dum%*%SGI%*%dum))/kk;   # f(x,y)
}}


## Draw the plot
par(mfrow=c(2,2));
image(xx,yy,zz,col=terrain.colors(20),xlab='X',ylab='Y',main='Terrain Colors');
segments(c(0,-4),c(-7,0),c(0,6),c(7,0),lty=2,col='black')
image(xx,yy,zz,col=heat.colors(20),xlab='X',ylab='Y',main='Heat Colors');
segments(c(0,-4),c(-7,0),c(0,6),c(7,0),lty=2,col='black')
image(xx,yy,zz,col=topo.colors(20),xlab='X',ylab='Y',main='Topo Colors'); 
segments(c(0,-4),c(-7,0),c(0,6),c(7,0),lty=2,col='black')
image(xx,yy,zz,col=rainbow(20),xlab='X',ylab='Y',main='Rainbow Colors'); 
contour(xx,yy,zz,nlevels=3,col='black',drawlabels=F,add=T)