Monte Carlo simulation on sphere: unbiased random steps

Question

Im doing a Metropolis Monte Carlo simulation with particles on a sphere and have a question concerning the random movement in a given time step.

I understand that to obtain a uniform distribution of random points on the sphere to begin with it is not enough to use the naive simplest way (use spherical coordinates with constant R and pick random angles theta and phi), but one has to for example use one of these methods: http://mathworld.wolfram.com/SpherePointPicking.html

Looking at another code for a Monte Carlo on a sphere I see a fairly complicated way to generate random moves: pick a random particle, calculate the rotation matrix moving it to the north pole, find a random cartesian vector less than a certain length and move it to the north pole, normalize the cartesian vector and then rotate it back to the vicinity of the original particle position.

This is all to get an unbiased new random position. I don`t understand the rationale completely although I suspect it is connected to detailed balance. But my feeling is that there should be an easier (i.e. faster) way to do this. Actually, intuitively I feel that in this case it is ok to find two small random angles theta and phi and add them to the position of the particle - or would this be a mistake?

Answer 1

The rationale behind the algorithm you describe is that the point is mapped to a pole, a tangent plane is generated, a Cartesian vector is generated uniformly in the usual way, this is then mapped back to the sphere, and the rotation is finally reversed. Since the tangent plane is a derivative, for large spheres and small steps this is a good approximation.

Can we do better? Or at least steal someone else's code? Maybe. Consider this: you can already generate a uniform distribution of random points on a sphere, thanks to that article on MathWorld.

Without loss of generality, you may consider your current point to be a pole of the sphere. If you use the MathWorld algorithm to generate a random point, you know that the resulting point has uniform probability of being at any azimuthal angle.

The problem is then reduced to finding a bearing on a sphere, given two points' phi/theta values and then generating a new point an incremental surface distance along that bearing's great circle route.

For this new version of the problem, you can get all the code you need by importing your favourite GIS or projection library and using a spherical model of the Earth (with appropriate scaling).

Many people seem to reference this site as having a ton of good math on spherical navigation.