####
#### Some example code for drawing an ellipse 
####


###
### Suppose we are dealing with the quadratic form
### (x -x0)'Q(x-x0) = k
###

my.ellipse <- function(Q,x0,k){
  v1 <- eigen(Q)$vec[,1]
  v2 <- eigen(Q)$vec[,2]
  e1 <- sqrt(1/eigen(Q)$val[1]*k)
  e2 <- sqrt(1/eigen(Q)$val[2]*k)
  theta <- (0:100)*2*pi/100
  x <- e1*v1[1]*cos(theta) + e2*v2[1]*sin(theta) + x0[1]
  y <- e1*v1[2]*cos(theta) + e2*v2[2]*sin(theta) + x0[2]
  cbind(x,y)
}


## plotting a circle

plot(my.ellipse(matrix(c(1,0,0,1),2,2),c(0,0),1),type="l")


## plotting an ellipse with axis lengths 3 and 1
plot(my.ellipse(matrix(c(1/9,0,0,1),2,2),c(0,0),1),type="l")