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

## Richard Gill has offered 10,000 Euroes to anyone who can simulate the N
## directions of angular momentum vectors appearing in equation (16) of this
## experimental proposal of mine: http://arxiv.org/abs/0806.3078. Here I am
## attempting to provide such N directions. They are given by the vectors 'u'
## 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=46.
## While this is by no means a perfect simulation of my model, it does meet
## all the stringent conditions set out by Richard Gill for his challenge.

## Since after the explosion the angular momentum vectors 'u' moving along
## the z direction will be confined to the x-y plane, a 2D simulation is good
## enough for my proposed experiment.

set.seed(9875)

## Measurement angles for directions 'a' and 'b' in the equatorial plane

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

K <- length(angles)

corrs <- numeric(K)  ## Container for correlations

L <- numeric(K)

N <- 10^4  ## Sample size. Next try 10^5, or even 10^6

z <- runif(N, -1, 1)

p <- sqrt(1 - z^2)

r <- runif(N, 0, pi)

s <- runif(N, 0, pi)

t <- runif(N, -1, +1)

w <- runif(N, -1, +1)

x <- p * cos(r/0.69) * sign(t)

y <- p * 1.21 * (-1 + (2/(sqrt(1 + (3 * (0.963 * s)/pi))))) * sign(w)

u <- rbind(x, y, z)  ## A 3xN matrix. The N columns of u represent the
## x, y, and z coordinates of points on an approximate 2-sphere, S^2

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'

        corrs[i] <- sum(sign(-ua) * sign(ub))/N

        # corrs[i] <- sum(sign(-ua))/N
    }
}

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 = 2)
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("black", "black"), 
    lty = 1, col = c("blue", "red"))

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'

        corrs[j] <- sum(sign(-ua) * sign(ub))/N

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

        L[j] <- N

    }
}

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 = 2)
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("black", "black"), 
    lty = 1, col = c("blue", "red"))

plot(angles * 180/pi, -cos(angles) - corrs, type = "l", col = "blue", ylim = c(-0.04, 
    +0.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(L), col = "green")
lines(angles * 180/pi, -1/sqrt(L), col = "green")