Surfaces.R

# Code modified from Zen.R

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))
Bell <- function(phi) {
    phi <- abs(phi)
    if (phi > pi) {
        phi <- 2 * pi - phi
    }
    return((phi - (pi/2))/(pi/2))
}
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")

This is the surface for a LHV model that is rotationally invariant (Bell surface):

BC = matrix(nrow = k, ncol = k, data = sapply(angles, function(t) sapply(t - 
    angles, Bell)), byrow = TRUE)
persp(x = angles, y = angles, z = BC, zlim = c(-1, 1), col = "green", theta = 135, 
    phi = 30, scale = FALSE, xlab = "alpha", ylab = "beta")

We now check overall performance by computing average distances, first between OM and Chistian's surface, then QM and Bell surface:

mean(abs(QM - corrs))

## [1] 0.1963

mean(abs(QM - BC))

## [1] 0.1349