ChaoticUnsharpBall3.R

richard — Feb 26, 2014, 10:50 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 45 degrees and 90 degrees give different interesting curves!
## Obvious fact: 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^4 (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^4


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

par(mfrow = c(3, 4))

for (R in c(45, 60, 70, seq(75, 87, by = 2), 88, 90)) {


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


    for (i in 2:(K-1)) {
        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
    }
    corrs[1] <- 1
    corrs[K] <- 0


    plot(angles * 180/pi, corrs, type = "l", col = "blue", 
        main = paste("R = ", R, " degrees"), 
        ylim = c(0, 1), xlim = c(0, 90), xlab = "", ylab = "",
        xaxp = c(0, 90, 3), yaxp = c(0, 1, 4))
    lines(angles * 180/pi, cos(angles), col = "black")      
}