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