## First Experiment
##
## Standard Bell experiment
##
## My code simulates 4000 trials of a Bell experiment in which:
##
## Alice uses random polarizer settings equal either to 0 deg or to 45 degrees
##
## Bob uses polarizer settings either equal to 22.5 degrees or to -22.5 degrees
##
## Each time the source shoots out a photon polarized at a random angle between 0 and 180 degrees
## towards Alice, and the same angle + 90 degrees towards Bob
##
## The detectors give +/-1 according to whether or not the cosine of twice the difference
## between the setting and the incoming polarisation is positive.
##
## I find three correlations equal to about -0.5 and one equal to about +0.5
set.seed(12345)
Avec <- sample(c(0, 45), replace = TRUE, size = 4000)
Bvec <- sample(c(22.5, -22.5), replace = TRUE, size = 4000)
LambdaVec <- runif(4000, 0 , 180)
Xvec <- rep(0, 4000)
Yvec <- rep(0, 4000)
for (i in 1:4000){
A <- Avec[i]
B <- Bvec[i]
ThetaA <- LambdaVec[i]
ThetaB <- ThetaA + 90
Xvec[i] <- ifelse((cos(2 * (A - ThetaA) * (2 * pi) / 360) > 0), +1, -1)
Yvec[i] <- ifelse((cos(2 * (B - ThetaB) * (2 * pi) / 360) > 0), +1, -1)
}
r11 <- mean((Xvec*Yvec)[Avec == 0 & Bvec == 22.5])
r12 <- mean((Xvec*Yvec)[Avec == 0 & Bvec ==-22.5])
r21 <- mean((Xvec*Yvec)[Avec == 45 & Bvec == 22.5])
r22 <- mean((Xvec*Yvec)[Avec == 45 & Bvec == -22.5])
r11; r12; r21; r22
## [1] -0.56
## [1] -0.5335893
## [1] -0.5262097
## [1] 0.484472
r22 - (r11 + r12 + r21)
## [1] 2.104271
## Second experiment
##
## Hoping to see the cosine curve?
##
## This time I keep Alice's setting fixed at 0 degrees.
##
## I let Bob's vary from 0 to 180.
## Each time I do 4000 trials so I can draw the whole "cosine curve"
## (which some people seem to expect) of the correlations,
## as the difference between the orientations varies from 0 to 180 degrees
set.seed(678910)
A <- 0
corr <- rep(0, 181)
i <- 1
for (B in 0:180){
ThetaA <- runif(4000, 0 , 180)
ThetaB <- ThetaA + 90
X <- ifelse((cos(2 * (A - ThetaA) * (2 * pi) / 360) > 0), +1, -1)
Y <- ifelse((cos(2 * (B - ThetaB) * (2 * pi) / 360) > 0), +1, -1)
corr[i] <- mean(X * Y)
i <- i + 1
}
plot(0:180, corr, type = "l")

## Third experiment
##
## This time I let Alice go from 0 to 90 degrees and hold Bob at 0
##
## Each time I do 4000 trials
set.seed(13243546)
B <- 0
corr <- rep(0, 91)
i <- 1
for (A in 0:90){
ThetaA <- runif(4000, 0 , 180)
ThetaB <- ThetaA + 90
X <- ifelse((cos(2 * (A - ThetaA) * (2 * pi) / 360) > 0), +1, -1)
Y <- ifelse((cos(2 * (B - ThetaB) * (2 * pi) / 360) > 0), +1, -1)
corr[i] <- mean(X * Y)
i <- i + 1
}
plot(0:90, corr, type = "l")
