I've been trying to generate a uniform spherical distribution in Python using uniform random sampling. For some reason my presumed spherical distribution looks more like an ovoid than like a sphere. I am using the fact that a sphere is defined: x^2 + y^2 + z^2 = R^2 and assuming R = 1 I get that the condition that the points should satisfy to be inside the sphere is: x^2 + y^2 + z^2 <= 1. For some reason this does not work. I have a perfect circular distribution from a top view (projection in the xy plane), but there is clearly a elliptical geometry in the planes xz and yz. ``` import numpy as np import numpy.random as rand import matplotlib.pyplot as plt
N = 10000
def sample_cube(number_of_trials):
"""
This function receives an integer and returns a uniform sample of points in a square.
The return object is a list of lists, which consists of three entries, each is a list
of the copordinates of a point in the space.
"""
cube = rand.uniform(low = -1, high = 1, size = (number_of_trials, 3))
x_cube = cube[:,0]
y_cube = cube[:,1]
z_cube = cube[:,2]
return [x_cube, y_cube, z_cube]
def sample_sphere(cube):
"""
This function takes a list on the form [x_cube, y_cube, z_cube] and then sample
an spherical distribution of points.
"""
in_sphere = np.where(np.sqrt(cube[0]**2 + cube[1]**2 + cube[2]**2) <= 1)
return in_sphere
"""
Main Code
"""
cube = sample_cube(N)
sphere = sample_sphere(cube)
print(sphere[0])
print(cube)
print(cube[0][sphere[0]])
x_in_sphere = cube[0][sphere[0]]
y_in_sphere = cube[1][sphere[0]]
z_in_sphere = cube[2][sphere[0]]
fig = plt.figure()
ax = plt.axes(projection = "3d")
ax.scatter(x_in_sphere, y_in_sphere, z_in_sphere, s = 1, color = "black")
plt.show()
plt.clf
![Distribution generated by the code above](https://i.stack.imgur.com/N83Ui.png)
![Side view of the distribution (plane xz)](https://i.stack.imgur.com/Di5kJ.png)
![Top view of distribution (plane xy)](https://i.stack.imgur.com/8uQvD.png)
I was simply trying to get a uniform sphere. There should be something wrong with the approach, but I can not spot the mistake.
Aucun commentaire:
Enregistrer un commentaire