My response to the 10,000 + 5,000 Euros challenge by Richard Gill:

## Richard Gill had offered 10,000 Euroes to anyone who can simulate the N
## directions of angular momentum vectors appearing in equation (16) of my
## experimental proposal: http://arxiv.org/abs/0806.3078. In this simulation
## I provided such N directions. They are given by the vectors 'v' in this
## simulation. He has also offered further 5,000 Euros to me if my proposed
## experiment is realized successfully. I am hopeful that that will happen
## some day. The details of these challanges by Richard Gill can be found
## here: http://www.sciphysicsforums.com/spfbb1/viewtopic.php?f=6&t=52#p1898.

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

## Since after the bomb explosion the angular momentum vectors 'v' moving
## along the z direction will be confined to the x-y plane, a 2D simulation
## is good enough for my proposed experiment: http://arxiv.org/abs/1211.0784.

set.seed(9875)

M <- 10^5  ## Size of the pre-ensemble. Next, try 10^6, or even 10^7

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

K <- length(angles)

corrs = matrix(nrow = K, ncol = K, data = 0)  ## Container for correlations

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)  ## A 3xM matrix. The M columns of u represent the x, y,
## and z coordinates of M uniformly distributed points on the 2-sphere, S^2.
## For a demonstration, see, for example, http://rpubs.com/gill1109/13340.

eta <- runif(M, 0, pi)

f <- -1 + (2/(sqrt(1 + (3 * eta/pi))))

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'

        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, j] <- sum(sign(va) * sign(-vb))/N

        ## corrs[j] <- sum(sign(vb))/N

    }
}

(N)

## [1] 66504

e <- cbind(o, p, q)

write.csv(e, file = "AliceDirections.csv")  ## Alice's N spin directions
write.csv(-e, file = "BobDirections.csv")  ## Bob's N spin directions

par(mar = c(0, 0, 2, 0))
persp(x = angles, y = angles, main = "Correlations predicted by my 3-sphere model agree exactly with the corresponding QM predictions", 
    z = corrs, zlim = c(-1, 1), col = "pink", theta = 135, phi = 30, scale = FALSE, 
    xlab = "alpha", ylab = "beta")

## We now calculate the four Bell-test correlations, with a and b fixed:

alpha <- 0 * pi/180
beta <- 45 * pi/180
a <- c(cos(alpha), sin(alpha), 0)
b <- c(cos(beta), sin(beta), 0)
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'
(E_0_45 <- sum(sign(va) * sign(-vb))/N)

## [1] -0.7079

(N)

## [1] 52236

alpha <- 0 * pi/180
beta <- 135 * pi/180
a <- c(cos(alpha), sin(alpha), 0)
b <- c(cos(beta), sin(beta), 0)
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'
(E_0_135 <- sum(sign(va) * sign(-vb))/N)

## [1] 0.7024

(N)

## [1] 52269

alpha <- 90 * pi/180
beta <- 45 * pi/180
a <- c(cos(alpha), sin(alpha), 0)
b <- c(cos(beta), sin(beta), 0)
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'
(E_90_45 <- sum(sign(va) * sign(-vb))/N)

## [1] -0.7062

(N)

## [1] 52281

alpha <- 90 * pi/180
beta <- 135 * pi/180
a <- c(cos(alpha), sin(alpha), 0)
b <- c(cos(beta), sin(beta), 0)
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'
(E_90_135 <- sum(sign(va) * sign(-vb))/N)

## [1] -0.7116

(N)

## [1] 52237

## The Bell-CHSH inequality is violated:

abs(E_0_45 - E_0_135 + E_90_45 + E_90_135)

## [1] 2.828

## It is worth noting that all four correlations calculated individually with
## fixed a and b lie on the correlation surface computed above with random a
## and b, thereby exhibiting the consistency of my model with the simulation.

## It is also worth noting that the number of initial states, N, in the last
## four calculations is the same within the margins of the square-root error:

(sqrt(N))

## [1] 228.6