ChaoticUnsharpBall2.R

richard — Feb 26, 2014, 8:53 AM

## Simulation of the Christian - Roth model 
## http://rpubs.com/chenopodium/joychristian

## ... and/or ...

## simulation of Caroline Thompson's chaotic spinning ball, various cap size
## http://freespace.virgin.net/ch.thompson1/Papers/The%20Record/TheRecord.htm

## Angular radius of the four circular caps is R

## Different values of R between pi/4 and pi/2 give different interesting curves!
## Conjecture: there is a convex combination of these curves which yields the cosine.

## This code is copy-paste from http://rpubs.com/chenopodium/joychristian
## My only changes:
## Decrease sample size M from 10^6 to 10^5 (for speed)
## Only look at angles between 0 and 90 degrees, 5 degree steps
## (save time, and increase resolution of plot: 
## magnification 2x vertically, 4x horizontally).


## Thanks especially to Joy and Chantal, but also to Michel, Daniel and others



set.seed(2938)


M <- 10^5

# R <- rbeta(M, 0.5, 0.5)
# s <- cos(R*pi/2)
# R <- R*180/pi

# s <- rbeta(M, 0.6, 0.4)
# R <- acos(s)
# R <- R*180/pi

gamma = 0.8
R <- pi/4 + gamma * pi/4
s <- cos(R)   ## Instead of acos((sin(Theta)^1.32)/3.16), Theta ~ Unif(0, pi/2)

R <- R * 180 / pi
R  ## radius in degrees (fixed instead of random)

[1] 81

# hist(R, breaks = 50,xlim=c(0, 100), freq = FALSE, 
#      main = "The cap radius distribution",
#      xlab = "Angular radius (degrees)")

angles <- seq(from = 0, to = 90, by = 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


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

e <- rbind(z, x)  ## 2 x M matrix


b <- c(cos(beta), sin(beta))  ## Measurement vector 'b'


for (i in 1:K) {
    alpha <- angles[i]
    a <- c(cos(alpha), sin(alpha))  ## Measurement vector 'a'
    ca <- colSums(e * a)  ## Inner products of cols of 'e' with 'a'
    cb <- colSums(e * b)  ## Inner products of cols of 'e' with 'b'
    good <- abs(ca) > s & abs(cb) > s  ## Select the 'states' 
    N <- sum(good)
    corrs[i] <- sum(sign(ca[good]) * sign(cb[good]))/N
    Ns[i] <- N
}


plot(angles * 180/pi, corrs, type = "l", col = "blue", main = "Two correlation functions", 
    xlab = "Angle (degrees)", ylab = "Correlation",
    ylim = c(0, 1))
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.2, legend = c(paste("Caroline, R = ", R, " degrees"), "cosine"), text.col = c("blue", 
    "black"), lty = 1, col = c("blue", "black"))