I've been trying to generate points on the surface of a sphere of radius "inner_radius", such that they're uniformly spread out. The algorithm works as expected for a radius of 1, but generates lesser than expected points for greater radii. I have looked through similar questions on here, but they seem to be for generating points throughout the volume and not just on the surface of the sphere.
import numpy as np
PI=np.pi
def spherical_to_cartesian(pol_ang,azim_ang,radius): #This function converts given spherical coordinates (theta, phi and radius) to cartesian coordinates.
return np.array((radius*np.sin(pol_ang) * np.cos(azim_ang),
radius*np.sin(pol_ang) * np.sin(azim_ang),
radius*np.cos(pol_ang))
)
def get_electron_coordinates_list(inner_radius,electron_count):
#Algorithm used was mostly taken from https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf . Explanations in code added by me.
electron_coordinate_list=[]
inner_area=4*(PI*inner_radius**2)
area_per_electron=inner_area/electron_count
pseudo_length_per_electron=np.sqrt(area_per_electron) #This is the side length of a square where the area of it is the area per electron on the sphere.
#Now, we need to get a value of angular space, such that angular space between electrons on latitude and longitude per electron is equal
#As a first step to obtaining this, we must make another value holding a whole number approximation of the ratio between PI and the pseudo_length. This will give the number of
#possible latitudes.
possible_count_of_lats=np.round(PI/pseudo_length_per_electron)
approx_length_per_electron_lat=PI/possible_count_of_lats #This is the length between electrons on a latitude
approx_length_per_electron_long=area_per_electron/approx_length_per_electron_lat #This is the length between electrons on a longitude
for electron_num_lat in range(int(possible_count_of_lats.item())): #The int(somenumpyvalue.item()) is used because Python cannot iterate over a numpy integer and it must be converted to normal int.
pol_ang=PI*(electron_num_lat+0.5)/possible_count_of_lats #The original algorithm recommended pol_ang=PI*(electron_num_lat+0.5)/possible_count_of_lats. The 0.5 appears to be added in order to get a larger number of coordinates.
#not sure if removing the 0.5 affects results. It didnt do so drastically, so what gives? Anyway, this gets the polar angle as PI*(latitudenumber)/totalnumberoflatitudes.
possible_count_of_longs=np.round(2*PI*np.sin(pol_ang)/approx_length_per_electron_long)
for electron_num_long in range(int(possible_count_of_longs.item())):
azim_ang=(2*PI)*(electron_num_long)/possible_count_of_longs #This gets the azimuthal angle as 2PI*longitudenumber/totalnumberoflongitudes
electron_coordinate=spherical_to_cartesian(pol_ang, azim_ang,inner_radius) #Converts the recieved spherical coordinates to cartesian so Manim can easily handle them.
electron_coordinate_list.append(electron_coordinate) #Add this coordinate to the electron_coordinate_list
print("Got coordinates: ",electron_coordinate) #Print the coordinate recieved.
print(len(electron_coordinate_list)," points generated.") #Print the amount of electrons will exist. Comment these two lines out if you don't need the data.
return electron_coordinate_list
get_electron_coordinates_list(1,100)
get_electron_coordinates_list(2,100)
Spherical_to_Cartesian() does nothing other than convert the spherical points to Cartesian.
For 100 points and radius 1, it generates 99 points. But, only 26 points are made if the radius is 2 and 100 points are requested.
If you can generate points uniformly in the sphere's volume, then to get a uniform distribution on the sphere's surface, you can simply normalize the vectors so their radius equals the sphere's radius.
Alternatively, you can use the fact that independent identically-distributed normal distributions are rotationally-invariant. If you sample from 3 normal distributions with mean 1 and standard deviation 0, and then likewise normalize the vector, it will be uniform on the sphere's surface. Here's an example:
import random
def sample_sphere_surface(radius=1):
x, y, z = (random.normalvariate(0, 1) for i in range(3))
scalar = radius / (x**2 + y**2 + z**2) ** 0.5
return (x * scalar, y * scalar, z * scalar)
To be absolutely foolproof, we can handle the astronomically unlikely case of a division-by-zero error when x, y and z all happen to be zero:
def sample_sphere_surface(radius=1):
while True:
try:
x, y, z = (random.normalvariate(0, 1) for i in range(3))
scalar = radius / (x**2 + y**2 + z**2) ** 0.5
return (x * scalar, y * scalar, z * scalar)
except ZeroDivisionError:
pass