fred.R

set.seed(1234)
M <- 1000000
theta <- runif(M, 0, 360)
beta1 <- 0.32
beta2 <- 0.938
phi <- 3
xi <- 0
s1 <- theta
s2 <- theta + 180
cosD <- function(x) cos(pi * x / 180)
sinD <- function(x) sin(pi * x / 180)
lambda1 <- beta1 * cos(theta/phi)^2
lambda2 <- beta2 * cos(theta/phi)^2

experiment <- function(a, b){
  Aa1 <- ifelse(abs(cosD(a - s1)) > lambda1, sign(cosD(a - s1)), 0)
  Aa2 <- ifelse(abs(cosD(a - s1)) < lambda2, sign(sinD(a - s1 + xi)), 0)
  Bb1 <- ifelse(abs(cosD(b - s2)) > lambda1, sign(cosD(b - s2)), 0)
  Bb2 <- ifelse(abs(cosD(b - s2)) < lambda2, sign(sinD(b - s2 + xi)), 0)
  products <- c(Aa1*Bb1, Aa2*Bb2)
  rate <- sum(abs(products))
  corr <- sum(products)/rate
  list(rate = rate, corr = corr)
}

rates <- rep(0, 181)
corrs <- rep(0, 181)
a <- 0
i <- 1

for (b in 0:180) {results <- experiment(a, b); rates[i] <- results$rate; corrs[i] <- results$corr; i <- i + 1}

plot(0:180, rates / (2 * M), type = "l", ylim = c(0, 1))

plot(0:180, corrs, type = "l")