Gisin2opt.R

richard — Feb 21, 2014, 1:36 PM

## Gisin and Gisin model (unsymmetrized version, for speed)
## http://arxiv.org/abs/quant-ph/9905018

## Generate 1 million particle pairs


set.seed(4567)

M <- 10^6

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

e <- rbind(x, y, z)

unifs <- runif(M)



angles <- seq(from = 0, to = 360, by = 7.5) * 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
b <- c(cos(beta), sin(beta), 0)  ## Measurement vector "b"


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

    ca <- colSums(e*a)
    cb <- colSums(e*b)

    good <- unifs < abs(ca)
    N <- sum(good)
    corrs[i] <- sum(sign(ca[good])*sign(cb[good]))/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 = 2)
lines(angles*180/pi, cos(angles), col = "black")
points(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2)

legend(x=0, y=-0.5, legend = c("Gisin^2", "cosine"), text.col=c("blue", "black"), lty=1, col = c("blue", "black"))