## Counting the number of accepted particle pairs in http://rpubs.com/jjc/84238


set.seed(1234)            ## for reproducibility
Angles = seq(from = 0, to = 360, by = 7.2) * 2 * pi/360
K = length(Angles)

Ns  = matrix(nrow = K, ncol = K, data = 0) # Container for no. events A and B

M = 10^4       ## Size of the pre-ensemble (for speed!)

r = runif(M, 0, 2*pi)
z = runif(M, -1, +1)
h = sqrt(1 - z^2)
x = h * cos(r)
y = h * sin(r)

e = rbind(x, y, z) 

s = runif(M, 0, pi)
f = -1 + (2/sqrt(1 + ((3 * s)/pi)))
g = function(u,v,s){ifelse(abs(colSums(u*v)) > f, colSums(u*v), 0)}
t = function(u,v,s){abs(sign(g(u,v,s))) > 0} 
n = function(Q,u1,v1,s1,u2,v2,s2){length(Q[t(u1,v1,s1) & t(u2,v2,s2)])}

for (i in 1:K) {
     alpha = Angles[i]
     a = c(cos(alpha), sin(alpha), 0)  # Measurement direction 'a'

   for (j in 1:K) {  
        beta = Angles[j]
        b = c(cos(beta), sin(beta), 0)  # Measurement direction 'b'
        A = +sign(g(a,e,s))  # Alice's measurement results A(a, e, s) = +/-1
        B = -sign(g(b,e,s))  # Bob's measurement results B(b, e, s) = -/+1
        P = length((A*B)[A != 0 & B != 0]) # The number of non-vanishing events observed 
        Ns[i,j] = P # Total number of non-vanishing +/- events observed by Alice and Bob
    }
}

par(mar = c(0, 0, 2, 0))
persp(x = Angles, y = Angles, expand = 4,
      main = "Relative frequency of simultaneous events (N/M), in [0, 1]", 
      z = Ns/M, zlim = c(0, 1), col = "bisque", theta = 135, phi = 30, 
      scale = FALSE, xlab = "alpha", ylab = "beta", zlab = "N / M")