Two methods sampling on a sphere:

Sample_On_Sphere.py

#!/usr/bin/python3
# -*- coding: utf-8 -*-
# 
# http://marc-b-reynolds.github.io/math/2018/06/21/SFPoints4ET.html

import numpy as np

"""
// constant turning rate: 
//   TX = cos(2pi K)
//   TY = sin(2pi K) 
//   K  = frac(phi) = 1/phi = (sqrt(5)-1)/2
const double TX = -0.73736887807831985597317725478205829858779907226562;
const double TY = -0.67549029426152362720614519275841303169727325439453;

typedef struct {
  double x,y;     // incrementally computed point on circle
  double z,dz;    // incrementally computed height on cap
} sf_walk_t;

// n = number of points to generate
// h = height of cap (ex: half-sphere=1, full-sphere=2)
void sf_walk_init(sf_walk_t* w, uint32_t n, float h)
{
  w->x  = 1.0;
  w->y  = 0.0;
  w->z  = 1.0;
  w->dz = h/n;
}

void sf_walk_next(sf_walk_t* w, float* v)
{
  double x=w->x, y=w->y;
  double ct,st;

  // current disc to cap mapping values
  ct   = w->z; 
  st   = sqrt(1-ct*ct);
  
  // output current point on cap
  v[0] = (float)(st*x);
  v[1] = (float)(st*y);
  v[2] = (float)(ct);
  
  // update point on circle: turn by 2pi*K
  w->x  = TX*x-TY*y;
  w->y  = TY*x+TX*y;

  // update height in cap position
  w->z -= w->dz;
}
"""

def Fib(n, r, h=2.0):
    """
    n = number of points to generate
    h = height of cap (ex: half-sphere=1, full-sphere=2)
    """
    TX = -0.73736887807831985597317725478205829858779907226562
    TY = -0.67549029426152362720614519275841303169727325439453
    wx, wy, wz, dz = 1.0, 0.0, 1.0, h/n

    pts = []
    for i in range(n):
        x = wx
        y = wy 
        ct = wz
        st = np.sqrt(1 - ct*ct)
        nx = st * x
        ny = st * y
        nz = ct
        wx = TX*x - TY*y
        wy = TY*x + TX*y
        wz -= dz 
        pts.append((nx, ny, nz))
    return pts 
    


# https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf

def Uniform(N, r = 1.0):
    points = []
    cnt = 0
    a = 4 * np.pi / N
    d = np.sqrt(a)
    M = int(np.round(np.pi / d))
    d_theta = np.pi / M
    d_phi = a / d
    for m in range(M):
        theta = np.pi * (m + 0.5) / M
        M_phi = int(round(2 * np.pi * np.sin(theta) / d_phi))
        for n in range(M_phi):
            phi = 2 * np.pi * n / M_phi
            x = r*np.sin(theta)*np.cos(phi)
            y = r*np.sin(theta)*np.sin(phi)
            z = r*np.cos(theta)
            points.append((x,y,z))
    return points

def Test_Uniform():
    import numpy as np
    import matplotlib.pyplot as plt    
    from mpl_toolkits import mplot3d

    pts = np.array(Uniform(100))
    ax = plt.axes(projection='3d')
    
    ax.scatter3D(pts[:, 0], pts[:, 1], pts[:, 2], cmap='Greens')
    plt.show()

    
def Test_Fib():
    import numpy as np
    import matplotlib.pyplot as plt    
    from mpl_toolkits import mplot3d
    pts = np.array(Fib(n = 1000, r = 10, h = 0.1))
    ax = plt.axes(projection='3d')    
    ax.scatter3D(pts[:, 0], pts[:, 1], pts[:, 2], cmap='Greens')
    plt.show()
    

if __name__ == '__main__':
    Test_Uniform()
    Test_Fib()