Uniformly distribute points in 4D

Question

I am trying to distribute points onto a 4D surface. I have an implementation that evenly distributes data points onto a 3D surface octant.

n_samples = 100
N = int(np.floor(np.sqrt(0.25+2*n_samples) -0.5))

phis = np.linspace(start=0, stop=0.25*np.pi, num=N)

# points = np.empty((int(N*(N+1)/2), 3))
points = []
for n in range(1,N+1):
    phi = phis[n-1]
    for theta in np.linspace(start=0, stop=0.5*np.pi, num=n):
        x = np.sin(phi)*np.cos(theta)
        y = np.sin(phi)*np.sin(theta)
        z = np.cos(phi)
        points.append([x,y,z])

points = np.array(points)

df = pd.DataFrame(points, columns=['x', 'y', 'z'])
fig = px.scatter_3d(df, x='x', y='y', z='z')
fig.show()

I am having difficulties in extending this to a 4D surface and visualizing it. Am I thinking about this correct? I need all my points evenly distributed in an octant. For 3D, I am leveraging the quadratic equation to distribute the points, for 4D, I am leveraging the cubic equation for the same.

n_samples = 100
poly = [1, 1.5, 0.5, 3*n_samples]
N = int(np.floor(-np.real(np.roots(poly)[0])) - 1)

phi_1s = np.linspace(start=0, stop=0.25*np.pi, num=N)

# points = np.empty((int(N*(N+1)/2), 3))
points = []
for n in range(1,N+1):
    phi_1 = phi_1s[n-1]
    for phi_2 in np.linspace(start=0, stop=0.25*np.pi, num=n):
        for phi_3 in np.linspace(start=0, stop=0.5*np.pi, num=n):
            x1 = np.cos(phi_1)
            x2 = np.sin(phi_1)*np.cos(phi_2)
            x3 = np.sin(phi_1)*np.sin(phi_2)*np.cos(phi_3)
            x4 = np.sin(phi_1)*np.sin(phi_2)*np.sin(phi_3)
            points.append([x1, x2, x3, x4])

points = np.array(points)

df = pd.DataFrame(points, columns=['x', 'y', 'z', 'c'])
fig = px.scatter_3d(df, x='x', y='y', z='z', color='c')

Any thoughts on this is appreciated..