pearle2.R

## A simulation of Pearle's (1970) model for the EPR-Bohm correlations
## In this script, I want to show how particle pairs move in and out of the
## sample as we increase one of the measurement angles, keeping the other
## fixed. We start with alpha = beta = 0 degrees, then increment beta by steps
## of 1 degree, till we come full circle at beta = 360 degrees.

set.seed(9875)       ## For reproducibility

## Uniform random points on sphere generated using 
## 'trig method' (method 3) of Dave Seaman:
## http://www.math.niu.edu/~rusin/known-math/96/sph.rand
## http://mathforum.org/kb/message.jspa?messageID=393612

## (a) Choose z uniformly distributed in [-1,1] 
## (b) Choose t uniformly distributed on [0, 2*pi)
## (c) Let r = sqrt(1-z^2) 
## (d) Let x = r * cos(t)
## (e) Let y = r * sin(t)

## Since the measurement directions will all be in the equatorial plane, 
## I will only generate z and x and treat them as x and y

## Measurement angles for setting 'a': directions in the equatorial plane

angles <- seq(from = 0, to = 360, by = 1) * 2 * pi/360

K <- length(angles)
corrs <- numeric(K)  ## Container for correlations
Ns <- numeric(K)  ## Container for number of states

beta <- 0 * 2 * pi/360  ## Measurement direction 'b' fixed, in equatorial plane

M <- 10^4

z <- runif(M, -1, 1)
t <- runif(M, 0, 2 * pi)
r <- sqrt(1 - z^2)
x <- r * cos(t)

e <- rbind(z, x)  ## 2 x M matrix
## The M columns of e represent the x and y coordinates
## of M uniform random points on the sphere S^2

U <- runif(M)
s <- (2/sqrt(3*U+1)) - 1  # Pearle's "r" is arc cosine of "s"; divided by pi/2

b <- c(cos(beta), sin(beta))  ## Measurement vector 'b'

## Loop through measurement vectors 'a'

good0 <- logical(M)

add <- integer(K)
leave <- integer(K)
inout <- matrix(TRUE, M, K)

for (i in 1:K) {
    alpha <- angles[i]
    a <- c(cos(alpha), sin(alpha))  ## Measurement vector 'a'
    ca <- colSums(e * a)  ## Inner products of cols of 'e' with 'a'
    cb <- colSums(e * b)  ## Inner products of cols of 'e' with 'b'
    good <- abs(ca) > s & abs(cb) > s  ## Select the 'states' 
    inout[ , i] <- good
    N <- sum(good)
    corrs[i] <- sum(sign(ca[good]) * sign(cb[good]))/N
    Ns[i] <- N
    add[i] <- sum(!good0 & good)
    leave[i] <- sum(good0 & !good)
    good0 <- good
}



plot(angles * 180/pi, corrs, type = "l", col = "blue", main = "Two correlation functions", 
    xlab = "Angle (degrees)", ylab = "Correlation")
lines(angles * 180/pi, cos(angles), col = "black")
legend(x = 0, y = -0.5, legend = c("Pearle", "cosine"), text.col = c("blue", 
    "black"), lty = 1, col = c("blue", "black"))

add[1] <- add[K]
leave[1] <- leave[K]

plot(angles * 180/pi, add, type = "l", main = "Numbers entering and leaving sample",
     sub = "Increment = 1 degree, total = 10^4",
     xlab = "Angle (degrees)",
     ylab = "Number", col = "blue")
lines(angles * 180/pi, leave, type = "l", col = "green")
legend(x = 45, y = 15, legend = c("Entering", "Leaving"), 
       text.col = c("blue", "green"), lty = 1, col = c("blue", "green"))

sum((add - leave)[-1])

## [1] 0

plot(angles * 180/pi, Ns, type = "l", col = "blue", 
     main = "Number of accepted particle pairs (out of 10^4)", 
     xlab = "Angle (degrees)", ylab = "Sample size",
     ylim = c(0, 10^4))

# inout <- as.logical(inout)
# dim(inout) <- c(M, K)

plot((0:360)[inout[1, ]], rep(1, 361)[inout[1, ]], pch =".", 
     xlim = c(0, 360), ylim = c(1,100),
     main = "100 particle pairs (plotted if in sample)",
     ylab = "index", xlab = "Angle (degrees)")
for (j in 2:100) points((0:360)[inout[j, ]], rep(j, 361)[inout[j, ]], pch =".", ylim = c(0,100))