richard — Feb 26, 2014, 3:53 PM
## Simulation of the Christian - Roth model
## http://rpubs.com/chenopodium/joychristian
## ... and/or ...
## simulation of Caroline Thompson's chaotic spinning ball, unsharp caps
## http://freespace.virgin.net/ch.thompson1/Papers/The%20Record/TheRecord.htm
## Angular radius of the four circular caps is a drawing from the distribution
## of R = (1 + U ^ gamma) * 45 degrees
## gamma = 0.46 is a little less than 1/2
## This particular choice obtained by least squares fit by RDG
## The code below: copy-paste from http://rpubs.com/chenopodium/joychristian
## Only changes:
## Increase sample size M from 10^6 to 10^7 (sqrt 10 higher precision)
## Only look at angles between 0 and 90 degrees
## (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(9082)
M <- 10^7
R <- (1 + runif(M)^(0.46)) * 45
hist(R, breaks=1000, xlim=c(45,90), freq=FALSE,
main = "The cap radius distribution",
xlab = "Angular radius (degrees)")
s <- cos(R * pi / 180)
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 N's
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")
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("Richard", "cosine"), text.col = c("blue",
"black"), lty = 1, col = c("blue", "black"))
angles[66] * 180 / pi
[1] 65
corrs[66]
[1] 0.4207
cos(angles[66])
[1] 0.4226
plot(angles * 180/pi, corrs - cos(angles), type = "l")
max(abs(corrs - cos(angles)))
[1] 0.002791
1/sqrt(Ns[66])
[1] 0.0004127