drawMahal <- 
function(x,x.mean,x.cov,whichlines=c(0.975,0.90,0.75),m=360,plot=TRUE,
         xlab="",ylab="",cex.lab=1.2, ...){
# draw ellipse corresponding to Mahalanobis distances
# PF, 06.12.2006
#
# x ... data
# x.mean ... mean of data
# x.cov ... cov-matrix of data
# whichlines: lines corresponding to quantiles of chi2
# m: number of data points on the ellipse
# plot ... if plot=TRUE the lines are plotted, otherwise not
#

# OUTPUT
# objects "mdX" and "mdY" with m rows and length(whichlines) columns
# (the columns refer to X- and Y-coordinates of one ellipse)

mdX=matrix(NA,nrow=m,ncol=length(whichlines))
mdY=matrix(NA,nrow=m,ncol=length(whichlines))

cov.svd <- svd(x.cov, nv = 0)
r <- cov.svd[["u"]] %*% diag(sqrt(cov.svd[["d"]]))

alphamd <- sqrt(qchisq(whichlines,2))

for (j in 1:length(whichlines)){
        e1md <- cos(c(0:(m-1))/m * 2 * pi) * alphamd[j]
        e2md <- sin(c(0:(m-1))/m * 2 * pi) * alphamd[j]
        emd <- cbind(e1md, e2md)
        ttmd <- t(r %*% t(emd)) + rep(1, m) %o% x.mean
	
	mdX[,j] <- ttmd[,1]
	mdY[,j] <- ttmd[,2]

	if (plot) lines(mdX[,j],mdY[,j], type = "l", ...)
}

list(mdX=mdX,mdY=mdY)
}