Dear Fred
Joy has already accepted the fact. You seem not to have noticed that he's changed the simulation model. He's also working with R, nowadays. I am not going to waste my time writing out a mathematical proof of something that Joy and I both agree on.
I suggest you install and learn R yourself. Both are pretty easy jobs.
Actually we are still only very close, still not yet spot-on. The problem is, what distribution to take for theta_0? At present Chantal Roth and Joy Christian together have homed in empirically on a pretty good choice, which in terms of Caroline Thompson's equivalent ball model, is circular caps of radius
R = acos((sin(Phi)^1.32)/3.16), Phi ~ Unif(0, pi/2).
Minkwe (lifted to 3-D - another of my innovations which Joy has taken on board) had
R = acos((sin(Phi)^2)/2), Phi ~ Unif(0, pi/2).
I've written an R script so one can experiment with various different probability distributions for R. And in particular, try various *constant* values of R. It is pretty clear that as you vary R between pi/4 and pi/2 you get various curves, none of which look much like a cosine, but such that a continuous convex linear combination of them can be found, which is exactly equal to the cosine! (OK: this is not a theorem, it is a conjecture, but it is a bloody good one, pardon my French!).
See:
http://rpubs.com/gill1109/ChaoticUnsharpBall2
The code is taken by copy-paste from Chantal and Joy's latest production. The only thing I have changed is the probability distribution of theta_0