> N <- 10^5
> x1 <- as.matrix(rnorm(N))
> x2 <- as.matrix(rnorm(N))
> x1 <- x1/sqrt(sum(x1^2))
> x2 <- x2/sqrt(sum(x2^2))
> (acos(sum(x1*x2))/pi)*180
[1] 89.94086

> library(MASS)
> rho <- 0.5
> Sig <-  matrix(c(1,rho,rho,1),2,2)
> X <- mvrnorm(N,mu=c(0,0),Sigma = Sig)
> x1 <- as.matrix(X[,1])
> x2 <- as.matrix(X[,2])
> x1 <- x1/sqrt(sum(x1^2))
> x2 <- x2/sqrt(sum(x2^2))
> (acos(sum(x1*x2))/pi)*180
[1] 59.90822

> rho <- -0.5
> Sig <-  matrix(c(1,rho,rho,1),2,2)
> X <- mvrnorm(N,mu=c(0,0),Sigma = Sig)
> x1 <- as.matrix(X[,1])
> x2 <- as.matrix(X[,2])
> x1 <- x1/sqrt(sum(x1^2))
> x2 <- x2/sqrt(sum(x2^2))
> (acos(sum(x1*x2))/pi)*180
[1] 119.8541