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]
  x <- ifelse(xunifs < cos(thetas - alpha)^2, 1, -1)
  y <- ifelse(yunifs < cos(thetas + pi/2 - beta)^2, 1, -1)
  corrs[i] <- mean(x * y)
}
plot(betas, corrs, type = "l", ylim = c(-1, +1))
abline(h = c(-1, 0, 1))