Random Rotation on a 3D sphere given an angle

Question

This question is in between computer graphic, probability, and programming, but since I am coding it for an Unity project in C# I decided to post it here. Sorry if not appropriate.

I need to solve this problem: given a object on a 3d sphere at a certain position, and given a range of degrees, sample points on the sphere uniformly within the given range.

For example: Left picture: the cube represents the center of the sphere, the green sphere is the starting position. I want to uniformly cover all surface of the circle within a certain degree, for example from -90 to 90 degrees around the green sphere. My approach (right picture) doesn't work as it over-samples points that are close to the starting position.

My sampler:

Vector3 getRandomEulerAngles(float min, float max)
{
    float degree = Random.Range(min, max);
    return degree * Vector3.Normalize(new Vector3(Random.Range(min, max), Random.Range(min, max), Random.Range(min, max)));
}

and for covering the top half of the sphere I would call getRandomEulerAngles(-90, 90).

Any idea?

Answer 1

Adapted from Nice Schertler, this is the code I am using

    Vector3 GetRandomAroundSphere(float angleA, float angleB, Vector3 aroundPosition)
    {
        Assert.IsTrue(angleA >= 0 && angleB >= 0 && angleA <= 180 && angleB <= 180, "Both angles should be[0, 180]");
        var v = Random.Range(0F, 1F);
        var a = Mathf.Cos(Mathf.Deg2Rad * angleA);
        var b = Mathf.Cos(Mathf.Deg2Rad * angleB);

        float azimuth = v * 2.0F * UnityEngine.Mathf.PI;
        float cosDistFromZenith = Random.Range(Mathf.Min(a, b), Mathf.Max(a, b));
        float sinDistFromZenith = UnityEngine.Mathf.Sqrt(1.0F - cosDistFromZenith * cosDistFromZenith);
        Vector3 pqr = new Vector3(UnityEngine.Mathf.Cos(azimuth) * sinDistFromZenith, UnityEngine.Mathf.Sin(azimuth) * sinDistFromZenith, cosDistFromZenith);
        Vector3 rAxis = aroundPosition; // Vector3.up when around zenith
        Vector3 pAxis = UnityEngine.Mathf.Abs(rAxis[0]) < 0.9 ? new Vector3(1F, 0F, 0F) : new Vector3(0F, 1F, 0F);
        Vector3 qAxis = Vector3.Normalize(Vector3.Cross(rAxis, pAxis));
        pAxis = Vector3.Cross(qAxis, rAxis);
        Vector3 position = pqr[0] * pAxis + pqr[1] * qAxis + pqr[2] * rAxis;
        return position;
    }