cs_

It would be interesting to view a skyplot of the GPS satellite trajectories, or at the least output the Dilution of Precision values comparing GPS with GPS+Galileo. This would help to illustrate the effect of additional satellites on measurement geometry, and thus on the precision of the output PVT solution.

A quick test with C++'s <random> library agrees. For n = 100,000,000 points, the rejection method beat out the improved spherical by 7.8%. Using the normal distribution was even slower.

  XYZ: 16.418057 seconds
  Gauss: 23.415560 seconds
  Spherical: 17.705553 seconds

I think that tends to produce a distribution clustered at the poles. At high latitudes, the impact of the longitude is much smaller than it would be around the equator. Mathworld has a nice write up that visualizes the effect and details more of the reasoning for the Gaussian method.

[0] http://mathworld.wolfram.com/SpherePointPicking.html