set.seed(12345)
alpha <- 0
gamma <- runif(10000, 0, 2 * pi)
results <- rep(0, 361)
for (i in 1:361){
beta <- (i - 1) * pi / 180
x <- 2 * ((gamma - alpha) %% (2 * pi) < pi / 2 | (alpha - gamma) %% (2 * pi) < pi / 2) - 1
y <- - (2 * ((gamma - beta) %% (2 * pi) < pi / 2 | (beta - gamma) %% (2 * pi) < pi / 2) - 1)
results[i] <- mean(x * y)
}
plot(0:360, results, type = "p", pch = 3, cex = 0.2,
main = "EPR correlations using a Gill-cosine rotor", xlab = "(a - b)", col = "blue")
points(0:360, cos(-pi/2 + results * pi/2), pch = 3, cex = 0.2, col = "red")
text(300, 0.7, "Bell's Boolean pairs", col = "blue")
text(300, 0.8, "Gill's cosine phase", col = "red")
