Plotting Points on a Sphere

Question

I am trying to generate random, uniformly distributed points on a sphere. However, the code is creating points that seem to create a disk instead. I believe that the issue resides in the "phirand" definition. Is the math there not correct? I used the same code in Matlab and it worked in that.

Code:

import numpy as np
import pylab
from scipy.integrate import odeint
import matplotlib.pyplot as plt
#import random
import mpl_toolkits.mplot3d.axes3d as p3
import random as rand

particlecount = 10 ## of particles to generate energies for energy generation
binsize = 15 #Determines bin size for historgram of electron energies
RStart = 0.02
phi1 = 0
phi2 = 180
phi1rad = phi1*(np.pi/180)
phi2rad = phi2*(np.pi/180)

#Generate random positions for each particle between s1 and s2
ICPositions = np.array([])
for i in range(0,particlecount):
    #In Spherical: Generates random position with boundaries of: S1<r<S2
    thetarand = (2*np.pi)*rand.uniform(0,1) #Random # generation for component y between s1 and s2
    phirand = np.arcsin((np.sin(phi2rad) - np.sin(phi1rad))*rand.uniform(0,1) + np.sin(phi1rad))
    xrand = RStart*np.sin(phirand)*np.cos(thetarand)
    yrand = RStart*np.sin(phirand)*np.sin(thetarand)
    zrand = RStart*np.cos(phirand)
    randArray = np.array([xrand,yrand,zrand])
    randArray = np.array(randArray,dtype = float)
    if ICPositions.size == 0:
        ICPositions = np.array([randArray])
    else:
        ICPositions = np.append(ICPositions,[randArray],axis = 0)

print(ICPositions)

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.scatter(ICPositions[:,0],ICPositions[:,1],ICPositions[:,2],c='r',marker='o')
ax.set_xlabel('x axis')
ax.set_ylabel('y axis')
ax.set_zlabel('z axis')
plt.show()

The `phirand` formula looks strange. It should be the same as `thetarand` (assuming you want uniform distribution over a shpere) — Marat, Feb 11 '20 at 02:46
Shouldn't be `phirand` defined like this? `phirand = np.arcsin(np.sin((phi2rad - phi1rad)*rand.uniform(0,1)) + np.sin(phi1rad))` — Ondro, Feb 11 '20 at 06:38
@Marat It is definitely different because theta has the range is [0,2pi] and phi range is [0,pi] — Tom, Feb 11 '20 at 07:57
@Ondro This does change the plot but the points follow a 3D parabolic shape instead of the disk. — Tom, Feb 11 '20 at 07:58
@Tom no, [0, pi] will result in only half the sphere. [-pi/2, pi/2] should work, but [0, 2pi] will work just the same — Marat, Feb 11 '20 at 12:34
@Marat If I set the range to be [-pi/2,pi/2] the points follow a 3D parabolic which is closer to what I want but still not it — Tom, Feb 12 '20 at 02:07
@Tom This question got addressed in one of your other [posts](https://stackoverflow.com/questions/60245192/half-of-matrix-size-missing/60293996#60293996). — Jenny, Feb 25 '20 at 22:03

Tom · Accepted Answer · 2020-02-27T23:32:39.443

I figured out the solution but do not understand why it works. Mathematically there is nothing wrong with my conversion definitions for x, y and z (After I Googled it and looked in a textbook). However, I saw another post where someone defined these coordinates slightly differently (Replacing sines with cosines for all phi components) without explaining why. The following worked when I replaced this:

xrand = RStart*np.sin(phirand)*np.cos(thetarand)
yrand = RStart*np.sin(phirand)*np.sin(thetarand)
zrand = RStart*np.cos(phirand)

With this:

xrand = RStart*np.cos(phirand)*np.cos(thetarand)
yrand = RStart*np.cos(phirand)*np.sin(thetarand)
zrand = RStart*np.sin(phirand)

Again I have no idea why this works but @Jenny gave a hint on another post that was asking a different question about the same code. Using [-90,90] instead of [0,180] is probably why this change is needed but again I am not certain why as sin[0,90,180] --> [0,1,0] while cos[-90,0,90] --> [0,1,0] and both cover the same numerical ranges. If someone has a solid mathematical/coding reason for this please comment on my answer to explain further.

Tom, you want a full sphere, or in other words, zrand has to run from [-1,1]. When defining phirand from [-90,90] the sine function will run from [-1,1] (but a cosine will only cover [0,1]). Similarly, if you define phirand = [0,180] the cosine will run from [-1,1]. — Jenny, Feb 28 '20 at 01:13
@Jenny I see now, thanks for the clarification. I think I made a few algebra mistakes. — Tom, Feb 28 '20 at 19:42

Plotting Points on a Sphere

1 Answers1