rmvtgauss <- function(n,d,location,scale,orientation)
              {matrix(rnorm(n*d),n,d) %*% 
                  diag(scale,d,d) %*% 
                  t(orientation) + 
                  matrix(location,n,d,byrow=TRUE)}

rmvtnorm <-function(n,mean,var) {rmvtgauss(n,
                                           length(mean),
                                           mean,
                                           sqrt(eigen(var)$values),
                                           eigen(var)$vectors)}