joy.R

## This is a tiny modification of Joy Christian's http://rpubs.com/jjc/16415
## which he claims is proof that he has won my challenge. He did not actually
## compute the four correlations which are needed to see if he has won or not,
## so here I do that for him, and then compute the CHSH quantity.
## Larger than 2.4, he wins; smaller than 2.4, I win.

######################################################################


## 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 setting 'a'' and 'b': directions in the equatorial
## plane

anglesAdeg <- c(0, 90)   ## CHSH angles for Alice, in degrees

anglesArad <- anglesAdeg * 2 * pi / 360   ## CHSH angles for Alice, in radians

anglesBdeg <- c(45, -45)  ## CHSH angles for Bob, in degrees

anglesBrad <- anglesBdeg * 2 * pi / 360

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

r <- runif(N, 0, pi)

s <- runif(N, 0, pi)

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

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

x <- cos(r/0.7) * sign(t)
y <- 1.2 * (-1 + (2/(sqrt(1 + (3 * (0.96 * s)/pi))))) * sign(w)

u <- rbind(x, y)  ## 2 x N matrix; N columns of u represent the
## x and y coordinates of points on an approximate circle.

corr <- matrix(0, 2, 2)

for (i in 1:2){
  alpha <- anglesArad[i]
  for (j in 1:2){
    beta <- anglesBrad[j]
    
    a <- c(cos(alpha), sin(alpha))  ## Measurement direction 'a'
    b <- c(cos(beta), sin(beta))  ## Measurement direction 'b'
    
    ua <- colSums(u * a)  ## Inner products of 'u' with 'a'
    ub <- colSums(u * b)  ## Inner products of 'u' with 'b'
    
    corr[i, j] <- sum(sign(-ua) * sign(ub))/N
    
    cat("\r", "alpha = ", anglesAdeg[i],
        ", beta = ", anglesBdeg[j], 
        ", correlation = ", corr[i, j], ", singletCorr = ", -cos(alpha - beta), sep = "")
    
  }
}

## 
alpha = 0, beta = 45, correlation = -0.7183, singletCorr = -0.7071068
alpha = 0, beta = -45, correlation = -0.72136, singletCorr = -0.7071068
alpha = 90, beta = 45, correlation = -0.2744, singletCorr = -0.7071068
alpha = 90, beta = -45, correlation = 0.28594, singletCorr = 0.7071068

CHSH <- corr[2, 2] - corr[1, 1] - corr[1, 2] - corr[2, 1]  ## typo corrected
CHSH

## [1] 2

cat("\r", "Christian wins is ", CHSH > 2.4, sep = "")

## 
Christian wins is FALSE