Zen.R

richard — May 4, 2014, 6:02 AM

angles = seq(from = 0, to = 360, by = 10) * pi / 180
k = length(angles)

corrs = matrix(nrow = k, ncol = k, data = 0)

N = 10^5
s = runif(N, 0, pi)
t = runif(N, 0, pi)
x = cos(s)
y = 1.2 * (-1 +  2 / (sqrt(1 + (3 * t / pi))))
e = rbind(x, y)

pb = txtProgressBar(min = 1, max = k * k, style = 3)
for (i in 1:k) {
  alpha = angles[i]
  a = c(cos(alpha), sin(alpha))
  for (j in 1:k) {
    beta = angles[j]
    b = c(cos(beta), sin(beta))
    corrs[i, j] = mean(sign(colSums(e * a)) * sign(colSums(-e * b)))
#    setTxtProgressBar(pb, (i - 1) * k + j)  ## Put back in if you run this live
  }
}

par(mfrow = c(2, 1))
par(mar = c(0, 0, 0, 0))
persp(x = angles, y = angles, z = corrs, zlim = c(-1, 1), col = "red", theta = 135, phi = 30, scale = FALSE, xlab = "alpha", ylab = "beta")
QM = matrix(nrow = k, ncol = k, data = sapply(angles, function(t) -cos(t - angles)), byrow = TRUE)
persp(x = angles, y = angles, z = QM, zlim = c(-1, 1), col = "blue", theta = 135, phi = 30, scale = FALSE, xlab = "alpha", ylab = "beta")