## 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")
