## A simulation of my 3-sphere model for the EPR-Bohm correlation. For a
## spectacular 2D-surface version of the simulaiton produced here, please see
## http://rpubs.com/jjc/16567. Further details and extensive discussion of
## the model can be found on my blog: http://libertesphilosophica.info/blog/.

## The theoretical description of the model can be found in this paper:
## http://arxiv.org/abs/1405.2355 (see also http://lccn.loc.gov/2013040705).

## This is S^2 version (aka 3D version: think of S^2 as a surface in R^3)
## This version has been adapted from Richard Gill's optimized version of
## Michel Fodje's original simulation of the model, which can be found here:
## http://rpubs.com/gill1109/EPRB3opt. Later Richard Gill improved his 3D
## version by employing the exact probability distribution derived by Philip
## Pearle in his classic 1970 paper: http://rpubs.com/gill1109/Pearle. It
## should be noted, however, that, unlike Pearle's model, the 3-sphere model
## has nothing whatsoever to do with data rejection or detection loophole.

## All of the above simulations are inspired by the original simulation of
## the 3-sphere model by Chantal Roth, https://github.com/chenopodium/JCS2.

set.seed(9875)

## For reproducibility. Replace integer seed by your own, or delete this line
## and let your computer dream up one for you (system time + process ID).

## 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 See, for example,
## http://rpubs.com/gill1109/13340 for an illustration in R.

## (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).

## 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^7  ## Size of the pre-ensemble. Next, try 10^8 for better results

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

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

eta <- runif(M, 0, pi)  ##  My initial eta_o, or Michel Fodje's 't'

f <- -1 + (2/sqrt(1 + ((3 * eta)/pi)))  ## Pearle's 'r' is arc cosine of 'f'

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

## Loop through measurement vectors 'a' (except last = 360 degrees = first)

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

    ua <- colSums(u * a)  ## Inner products of 'u' with 'a'
    ub <- colSums(u * b)  ## Inner products of 'u' with 'b'

    good <- abs(ua) > f & abs(ub) > f  ## Sets the topology to that of S^3

    o <- x[good]
    p <- y[good]
    q <- z[good]

    v <- rbind(o, p, q)  ## N spin directions produced at the source

    g <- f[good]  ## The initial state of spin is defined by the pair (v, g)

    N <- length(v)/3  ## The number of complete states (v, g) at the source

    va <- colSums(v * a)  ## Inner products of 'v' with 'a'
    vb <- colSums(v * b)  ## Inner products of 'v' with 'b'

    corrs[i] <- sum(sign(va) * sign(-vb))/N

    Ns[i] <- N
}

corrs[K] <- corrs[1]
Ns[K] <- Ns[1]

plot(angles * 180/pi, corrs, type = "l", col = "blue", main = "Two correlation functions", 
    xlab = "Angle (degrees)", ylab = "Correlation")
points(angles * 180/pi, corrs, col = "blue", pch = ".", cex = 4)
lines(angles * 180/pi, corrs, col = "blue")
lines(angles * 180/pi, -cos(angles), col = "red")
points(angles * 180/pi, -cos(angles), col = "red", pch = ".", cex = 1)

legend(x = 0, y = 0.9, legend = c("S^3", "-cos"), text.col = c("blue", "red"), 
    lty = 1, col = c("blue", "red"))

plot(angles * 180/pi, -cos(angles) - corrs, type = "l", col = "blue", ylim = c(-5e-04, 
    +5e-04), main = "Difference between the -cos (red) and S^3 (blue) correlations, with +/-1 standard error in green")
abline(h = 0, col = "red", lwd = 1)
lines(angles * 180/pi, 1/sqrt(Ns), col = "green")
lines(angles * 180/pi, -1/sqrt(Ns), col = "green")