# Ways to generate random points on a sphere¶

In [1]:
#All necesary import
import sys
import time

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline

sys.path.append('../')

from ten.utils import generate_random_points_in_sphere


## On the surface¶

The intuitive way but... don't work

Reference :

In [3]:
n = 3000
R = 5

theta = np.random.uniform(low = 0, high = 2*np.pi, size = n)
phi = np.random.uniform(low = 0, high = np.pi, size = n)

x = np.sin(phi)*np.cos(theta)
y = np.sin(phi)*np.sin(theta)
z = np.cos(phi)

fig = plt.figure(figsize=(12,12))
ax.scatter(x, y, z)

Out[3]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa775223ba8>

Note the large number of dots in the poles, there is no uniform distriuciÃ³n.

### The right way¶

In [5]:
n = 3000
theta = 2*np.pi*np.random.uniform(size=n)
phi = np.arccos(2*np.random.uniform(size=n) - 1)

x = np.sin(phi)*np.cos(theta)
y = np.sin(phi)*np.sin(theta)
z = np.cos(phi)

fig = plt.figure(figsize=(12,12))
ax.scatter(x, y, z)

Out[5]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa7750cf4a8>

Now we have a uniform distribution.But, what if we want to generate inside the sphere?

## Inside the sphere¶

### The wrong way¶

In [6]:
n = 3000
theta = 2*np.pi*np.random.uniform(size=n)
phi = np.arccos(2*np.random.uniform(size=n) - 1)

#theta = np.random.uniform(low=0, high=2*np.pi, size=n)
#phi = np.arcsin(2 * np.random.uniform(low=0, high=2*np.pi, size=n) - 1)
u = np.random.uniform(low=0, high=10, size=n)

x = np.sin(phi)*np.cos(theta)*u
y = np.sin(phi)*np.sin(theta)*u
z = np.cos(phi)*u

fig = plt.figure(figsize=(12,12))
ax.scatter(x, y, z)

Out[6]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa77510c4a8>

Obviously there is not a uniform distribution

### The trivial way, with a if...¶

In [7]:
n = 3000
R = 5
a = 0

x_l = []
y_l = []
z_l = []

t = time.time()

while a <= n:
x = np.random.uniform(low=-R, high=R, size=1)
y = np.random.uniform(low=-R, high=R, size=1)
z = np.random.uniform(low=-R, high=R, size=1)
if x*x + y*y + z*z <= R*R:
a += 1
x_l.append(x)
y_l.append(y)
z_l.append(z)

print("time =",time.time() - t)

fig = plt.figure(figsize=(12,12))

time = 0.30646824836730957

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa774f50780>