thetas <- runif(100000, 0, pi)
xunifs <- runif(100000)
yunifs <- runif(100000)
alpha <- 0
betas <- seq(from = 0, to = pi, length = 1000)
corrs <- rep(0, 1000)
for(i in 1:1000){
  beta <- betas[i]
  xprobs <- abs(cos(thetas - alpha))
  x <- 2  * (xunifs < xprobs) - 1
  yprobs <- abs(cos(thetas - beta))
  y <- 2  *  (yunifs < yprobs) - 1
  corrs[i] <- mean(x * y)
}
plot(betas, corrs, type = "l", ylim = c(-1, +1))
abline(h = c(-1, 0, +1))