Evenly distributed normal distribution on surface of a unit sphere Hello I am using Mathematica code to to generate NPoints of randomly generated points normally distributed around one point on the surface of a sphere using 
Return[Table[
   RandomVariate[NormalDistribution[phi, width]], {k, 1, NPoints}]];
to find my phi and
Return[Table[
   RandomVariate[NormalDistribution[theta, width]], {j, 1, NPoints}]];
to find my theta.
There is an issue when I plot the points though. They seem to smear along some direction when distributed around certain thetas and phis. I think the issue is the I need a correction term for projecting my points from a plane onto a sphere but I'm having a lot of trouble finding the term because I find changing coordinate systems difficult.
Any help would be greatly appreciated!
 A: I get the impression that perhaps you don't have a precise definition of "points normally distributed around one point on the surface of a sphere" in mind and are just looking for a pragmatic solution that looks like points normally distributed in a plane for deviations small compared to the radius of the sphere. If so, you might want to consider choosing the geodesic distance according to the distribution of the distance from the origin of points normally distributed in the plane. The probability density function for this distribution is proportional to $r\mathrm e^{-\beta r^2}$, where $r$ is the distance from the origin in the plane, so the cumulative distribution function is proportional to $\mathrm e^{-\beta r^2}$, so you can generate a uniform random variate $u$ and take
$$r = \sqrt{-\frac{\log u}\beta}$$
for the geodesic distance. Since the distance is limited to $[0,\pi]$, $u$ has to be uniformly distributed on $[\mathrm e^{-\beta\pi^2},1]$.
To get points distributed around the north pole, choose $\theta=r$ and $\phi$ uniformly in $[0,2\pi]$; then to get points distributed around some other point, just apply a suitable rotation.
This is computationally rather simple; however, it only works well for deviations small compared to the radius. If the probability is appreciable for geodesic distances greater than $\pi/2$, you'll get bunching towards the point's antipode, since the density of points at a certain geodesic distance from the point starts to decrease at $\pi/2$ and becomes zero again at the antipode.
Alternatively, you could choose points distributed according to the uniform distribution on the sphere weighted by $\mathrm e^{-\beta r^2}$, where $r$ is the geodesic distance from the point. That corresponds to replacing the factor $r$ by $\sin r$, so you'd need to sample from a distribution with probability density function proportional to $\sin r\mathrm\,e^{-\beta r^2}$. Integrating that leads to the imaginary error function. If you want to avoid complications from that, you could use rejection sampling, sampling from the above "uncorrected" distribution and accepting the samples with probability $\operatorname{sinc} r=\sin r/r$. For small to moderate deviations, this would be quite efficient, but even in the limit $\beta\to0$ the acceptance rate would still be
$$\frac{\displaystyle\int_0^\pi\frac{\sin r}r r\,\mathrm dr}{\displaystyle\int_0^\pi r\,\mathrm dr}=\frac2{\pi^2/2}=\left(\frac2\pi\right)^2\approx0.4\;.$$
