<p>Random points on the surface of a sphere can be generated by allowing the azimuthal angle $\phi$ to take a uniformly distributed random value between 0 and $2\pi$. To convert this to RA in degrees you multiply by $180/\pi$. To convert to hours, minutes and seconds you divide the $\phi$ in degrees by 15, which gives the hours, divide the remainder by 60 which gives the minutes and then divide the remainder of that by 60 to give the seconds. </p>

<p>The declination angle $\delta$ is not uniformly distributed, because the area covered at an angle $\delta$ goes as $\cos \delta$. </p>

<p>Note the potential confusion here between declination and the polar angle $\theta$, which is conventionally measured <em>down</em> from the pole. The area of a strip of width $\Delta \theta$ at a given $\theta$ on a unit sphere is $$\Delta A = 2\pi  \sin \theta\ \Delta \theta = 2\pi  \cos (\pi/2 - \delta)\ \Delta \delta = 2\pi  \cos \delta\ \Delta \delta$$</p>

<p>To translate this into something we can use, we say that the probability of finding a star in a small range of declinations
$$ dP \propto \cos \delta\ d\delta,$$
so the cumulative probability of finding a value of $\delta$ between $-90^{\circ}$ and $\delta$ is
$$ P(\delta) = \frac{\int^{\delta}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}}{\int^{+90^{\circ}}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}},$$
where the denominator makes sure the proportionality is correctly normalised. From this, we obtain
$$P(\delta) = \frac{1}{2}(1 + \sin \delta) \tag*{(1)}$$</p>

<p>The procedure for then obtaining a random declination is to assign a uniform random number to the cumulative probability $P$ between 0 and 1. Then from equation (1) we have
$$ \sin \delta = 2P -1$$
$$\delta = \sin^{-1}(2P-1),$$
where $P$ is a random number between 0 and 1.</p>

CC BY-SA 3.0
2017-08-29T14:57:50.787Z
@barrycarter Now sorted out - see edit. My original statement was correct.
ProfRob
CC BY-SA 4.0
2021-07-10T16:25:19.930Z
I've just done a Sage / Numpy version that you can rotate, pan & zoom. I'll post the link in the next comment. (Sorry that the code's a bit cryptic, I had to trim it down to make it small enough to fit).
PM 2Ring
CC BY-SA 4.0
2021-07-10T16:25:31.847Z
[Stars](https://sagecell.sagemath.org/?z=eJxtUstu2zAQvPMrFukhpMvYstwARQABOSXXAOndYCXaJiqRBB9una_vrOSiLtCLsDua2R3tyE0xpEK-TvFCJpOPIvkjdSjWyfghTOvBHkwdyx64dL7IL61SIjpwolt7qcRQTmi2j7QCQhvafm2EgIp-XT5kMpoG26snQZRA60OW3HNrS02e0oqxZJRGmZ1fSi5mnqiazhCeTZL3VZ_vlaiwBRBP4BLvG43NqxYyedb0EN2mZWTTwmhkoyaZyZbk-n0cQ9kN4p3RKNkhz4fy71BNIZrelUvXrLdKvIKa6ySvdP4g5jNT0SEkSoacp4xzHa28WsFHl5Ni8ed_1HW5xrJukaO_0UN24__PGPGMy9tk-jIfdi9zMSl3jw3WfTfj2H1L1Wo6JjfM5XJuM2e0ahEMwrvmuUj5_ry4Q-i49H8p9EBb5r3xmAADu0Hit-jDWCe_B6f_IW8jnlP7sN1OI-UxpO7uU9M0dzzCHWabT5iFe7wvCLu9Iq9A8in8lG-aDpxV92LGbNVvRxXJQQ==&lang=sage)
PM 2Ring
CC BY-SA 4.0
2022-04-06T23:39:11.623Z
@ProfRob Bravo! Thanks for taking the time towork this out so explicitly. Spherical coordinates tend to start with $\theta = 0$ at the top while latitude and declination put it in the middle, so I always get confused and just start trying $\sin^{-1}$ and $\cos^{-1}$ until something works.
uhoh
CC BY-SA 3.0
2017-08-29T11:09:52.617Z
2017-08-29T15:50:44.580Z
2017-08-29T15:50:44.580Z
ProfRob
<p>Supplemental answer for future readers here who are more interested in the &quot;uniformness&quot; than the &quot;randomness&quot; or &quot;realism&quot; of the distribution.</p>
<p>Here are two Stack Overflow answers that are similar, but the offset in the sunflower method can be tweaked to optimize the &quot;uniformness&quot; of the nearly uniform distribution.</p>
<ul>
<li><a href="https://stackoverflow.com/a/44164075/3904031">https://stackoverflow.com/a/44164075/3904031</a></li>
<li><a href="https://stackoverflow.com/a/26127012/3904031">https://stackoverflow.com/a/26127012/3904031</a></li>
</ul>
<p>The main reason that the patterns look different is that the &quot;poles&quot; where the spiral begins and ends is on the side for the &quot;fibonacci_sphere&quot; and the top for the &quot;sunflower_on_sphere&quot;.</p>
<p><a href="https://i.sstatic.net/oTctc.jpg" rel="noreferrer"><img src="https://i.sstatic.net/oTctc.jpg" alt="two nearly uniform distributions of points on a sphere" /></a></p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def sunflower_on_sphere(n=1000):
    # https://stackoverflow.com/a/44164075/3904031

    indices = np.arange(n) + 0.5

    phi = np.arccos(1 - 2 * indices / n)
    theta = np.pi * (1 + np.sqrt(5)) * indices

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

    return np.vstack([x, y, z])

def fibonacci_sphere(n=1000):
    # https://stackoverflow.com/a/26127012/3904031
    
    phi = np.pi * (3. - np.sqrt(5.))  # golden angle in radians

    theta = phi * np.arange(n)  # golden angle increment

    y = np.linspace(1, -1, n)  # y goes from 1 to -1
    
    r = np.sqrt(1 - y**2)  # radius at y
    
    x, z = [r*f(theta) for f in (np.cos, np.sin)]

    points = np.vstack([x, y, z])

    return points

n = 2000

xf, yf, zf = fibonacci_sphere(n)
xs, ys, zs = sunflower_on_sphere(n)

fig = plt.figure(figsize=[14, 7.5])

axf  = fig.add_subplot(1, 2, 1, projection='3d', proj_type = 'ortho')
axs  = fig.add_subplot(1, 2, 2, projection='3d', proj_type = 'ortho')

axf.scatter(xf, yf, zf, marker='.', c=yf, cmap='gray')
axs.scatter(xs, ys, zs, marker='.', c=xs, cmap='gray')

for ax in (axf, axs):
    ax.set_xlim(-1.1, 1.1)
    ax.set_ylim(-1.1, 1.1)
    ax.set_zlim(-1.1, 1.1)
    ax.set_box_aspect([1,1,1])
    ax.set_box_aspect([1,1,1]) # https://stackoverflow.com/a/68242226/3904031

axf.view_init(6, -60)
axs.view_init(6, -150)
axf.set_title('fibonacci_sphere(n='+str(n)+')')
axs.set_title('sunflower_on_sphere(n='+str(n)+')')
plt.show()
</code></pre>

CC BY-SA 4.0
2022-03-22T23:28:06.057Z
2022-03-22T23:28:06.057Z
uhoh
<p>This is really more computer science or mathematics than astronomy, but here's a function <code>random_point_on_unit_sphere()</code>  that creates random points on a unit 3D sphere, which is then used by another function <code>random_point_on_sky()</code> to convert to RA and dec values, using <a href="http://www.astropy.org" rel="nofollow noreferrer">astropy</a>. Finally, the function <code>print_random_star_coords()</code> print a number of RA,dec pairs.</p>

<pre><code>import numpy as np
from astropy.coordinates import SkyCoord
from astropy import units as u

def random_point_on_unit_sphere():
    while True:
        R   = np.random.rand(3) #Random point in box
        R   = 2*R - 1
        rsq = sum(R**2)
        if rsq &lt; 1: break       #Use r only if |r|&lt;1 in order not to favor corners of box
    return R / np.sqrt(rsq)     #Normalize to unit vector

def random_point_on_sky():
    p     = random_point_on_unit_sphere()
    r     = np.linalg.norm(p)
    theta = 90 - (np.arccos(p[2] / r)    / np.pi * 180)            #theta and phi values in degrees
    phi   =       np.arctan(p[1] / p[0]) / np.pi * 180
    c     = SkyCoord(ra=phi, dec=theta, unit=(u.degree, u.degree)) #Create coordinate
    return c.ra.hms, c.dec.deg                                     #Many different formats are possible, e.g c.ra.hour for decimal hour values

def print_random_star_coords(nstars):
    for n in range(nstars):
        RA,dec = random_point_on_sky()
        print(RA,dec)
</code></pre>

CC BY-SA 3.0
2017-08-29T10:15:15.670Z
2017-08-29T11:27:25.407Z
2017-08-29T11:27:25.407Z
pela
<p>Heres more python than you can shake a telescope at. I just used <a href="https://i.sstatic.net/G1rdk.png" rel="nofollow noreferrer">@ProfRob's algorithm</a>. This is just a python script, the real answer to the question is <a href="https://i.sstatic.net/G1rdk.png" rel="nofollow noreferrer">@ProfRob's answer</a> and I've just scripted it. The mathematics behind generating statistically uniform distributions is explained very nicely there as well.</p>
<p>Python is below the plots. You can see on an X-Y graph they thin out at the top and bottom.</p>
<p>The 3D plot is live, you can use your cursor and move the sphere around. Of course not HERE :-), but in your own python window if you run the script.</p>
<p><a href="https://i.sstatic.net/G1rdk.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/G1rdk.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/avktz.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/avktz.png" alt="enter image description here" /></a></p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]

degs, rads = 180/pi, pi/180

# do radians first, then convert later

nstars = 2000

ran1, ran2 = np.random.random(2*nstars).reshape(2, -1)

RA  = twopi * (ran1 - 0.5)
dec = np.arcsin(2.*(ran2-0.5)) # Hey Barry Carter!

funcs = [np.cos, np.sin]

cosRA,  sinRA  = [f(RA)  for f in funcs]
cosdec, sindec = [f(dec) for f in funcs]

x = cosRA * cosdec
y = sinRA * cosdec
z = sindec

decs = (np.arange(11)-5) * halfpi / 6.
RAs  = (np.arange(12)-5) * halfpi / 6.

theta = np.linspace(0, twopi, 101)

costh, sinth = [f(theta) for f in funcs]
zerth = np.zeros_like(theta)

fig = plt.figure()

ax  = fig.add_subplot(1, 1, 1, projection='3d')

ax.plot(x, y, z, '.k')

# lines of declination
xvals = costh * np.cos(decs)[:, None]
yvals = sinth * np.cos(decs)[:, None]
zvals = zerth + np.sin(decs)[:, None]
for x, y, z in zip(xvals, yvals, zvals):
    plt.plot(x, y, z, '-g', linewidth=0.8)

# lines of Right Ascention
xvals = costh * np.cos(RAs)[:, None]
yvals = costh * np.sin(RAs)[:, None]
zvals = sinth + np.zeros_like(RAs)[:, None]
for x, y, z in zip(xvals, yvals, zvals):
    plt.plot(x, y, z, '-r', linewidth=0.8)

ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-1.1, 1.1)
ax.set_zlim(-1.1, 1.1)

ax.view_init(elev=30, azim=15)

plt.show()

if True:
    plt.figure()
    plt.plot(degs*RA, degs*dec, '.k')
    plt.xlim(-180, 180)
    plt.ylim(-90, 90)
    plt.show()
</code></pre>

CC BY-SA 3.0
2017-08-29T13:23:13.100Z
Thanks a lot, so you are sure about the dec = np.arcsin(2.*(ran2-0.5)) ? It seems to be good but I don't understand very well why arcsin and not arcos for instance.
Richard
CC BY-SA 3.0
2017-08-29T13:53:22.507Z
Very nice! But when I try to rotate it only moves a little bit in the direction of my movement, after which I have to click/drag again for another little bit of movement. Aren't you supposed to be able to drag the plot around smoothly? Is it just me?
pela
CC BY-SA 3.0
2017-08-29T15:35:39.447Z
@pela it works nicely in my IDLE installation which I love, but the 3D plot response in my anaconda installation does not work well at all. This would be a good question for stackoverflow! You'd just edit the script down to the [Minimal, Complete, and Verifiable example](https://stackoverflow.com/help/mcve) before posting. You could link here also, but make sure to get rid of as much as possible.
uhoh
CC BY-SA 3.0
2017-08-29T18:06:38.787Z
@uhoh: Interesting, I have Anaconda install as well. Good idea with the Stackoverflow, thanks!
pela
CC BY-SA 4.0
2017-08-29T12:03:23.933Z
2022-04-06T23:39:26.363Z
2022-04-06T23:39:26.363Z
uhoh
<p>How to generate a uniform distribution of stars on the sky ? What I want is make a simulation of random points following uniform distribution on a part of the sky (making the assuption that we have the same amount of stars in any direction). The problem is how I do that because of the cos(dec) on the RA ?</p>
<p>What I did in python is this</p>
<pre><code>import numpy as np
RA = np.random.uniform(-180,180)
Dec = np.random.uniform(-90,90)
</code></pre>
<p>But it is totally false, since we will have more points on the poles (if I do many points, this will clearly not be uniform)...</p>
<p>Thanks a lot in advance!</p>

CC BY-SA 3.0
2017-08-29T09:08:29.767Z
What about `DEC=Arccos(np.random.uniform(-90,90))`?
J. Chomel
CC BY-SA 3.0
2017-08-29T09:44:04.447Z
Thank for the edit ! No I don't think so because arccos is defined between -1 and 1, but it is probably something in this direction.
Richard
CC BY-SA 3.0
2017-08-29T10:10:48.240Z
Something like dec=np.arcsin(np.random.uniform(-1,1)) looks a bit better, but is there someone to that can confirm me if it is mathematically correct or not ?
Richard
<positional-astronomy><python><mathematics>
Generate a uniform distribution on the sky
Richard
# Generate a uniform distribution on the sky

(Asked by: username_0 on 2017-08-29)

How to generate a uniform distribution of stars on the sky ? What I want is make a simulation of random points following uniform distribution on a part of the sky (making the assuption that we have the same amount of stars in any direction). The problem is how I do that because of the cos(dec) on the RA ?
What I did in python is this
import numpy as np
RA = np.random.uniform(-180,180)
Dec = np.random.uniform(-90,90)

But it is totally false, since we will have more points on the poles (if I do many points, this will clearly not be uniform)...
Thanks a lot in advance!

**Tags:** positional-astronomy, python, mathematics
> **Comments:**
> **username_1**: What about `DEC=Arccos(np.random.uniform(-90,90))`?  
> **username_0**: Thank for the edit ! No I don't think so because arccos is defined between -1 and 1, but it is probably something in this direction.  
> **username_0**: Something like dec=np.arcsin(np.random.uniform(-1,1)) looks a bit better, but is there someone to that can confirm me if it is mathematically correct or not ?  

---

## 4 answers

### Answer by username_4 on 2017-08-29. Score: 6 | Accepted answer

Heres more python than you can shake a telescope at. I just used @ProfRob's algorithm. This is just a python script, the real answer to the question is @ProfRob's answer and I've just scripted it. The mathematics behind generating statistically uniform distributions is explained very nicely there as well.
Python is below the plots. You can see on an X-Y graph they thin out at the top and bottom.
The 3D plot is live, you can use your cursor and move the sphere around. Of course not HERE :-), but in your own python window if you run the script.


import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]

degs, rads = 180/pi, pi/180

# do radians first, then convert later

nstars = 2000

ran1, ran2 = np.random.random(2*nstars).reshape(2, -1)

RA  = twopi * (ran1 - 0.5)
dec = np.arcsin(2.*(ran2-0.5)) # Hey Barry Carter!

funcs = [np.cos, np.sin]

cosRA,  sinRA  = [f(RA)  for f in funcs]
cosdec, sindec = [f(dec) for f in funcs]

x = cosRA * cosdec
y = sinRA * cosdec
z = sindec

decs = (np.arange(11)-5) * halfpi / 6.
RAs  = (np.arange(12)-5) * halfpi / 6.

theta = np.linspace(0, twopi, 101)

costh, sinth = [f(theta) for f in funcs]
zerth = np.zeros_like(theta)

fig = plt.figure()

ax  = fig.add_subplot(1, 1, 1, projection='3d')

ax.plot(x, y, z, '.k')

# lines of declination
xvals = costh * np.cos(decs)[:, None]
yvals = sinth * np.cos(decs)[:, None]
zvals = zerth + np.sin(decs)[:, None]
for x, y, z in zip(xvals, yvals, zvals):
    plt.plot(x, y, z, '-g', linewidth=0.8)

# lines of Right Ascention
xvals = costh * np.cos(RAs)[:, None]
yvals = costh * np.sin(RAs)[:, None]
zvals = sinth + np.zeros_like(RAs)[:, None]
for x, y, z in zip(xvals, yvals, zvals):
    plt.plot(x, y, z, '-r', linewidth=0.8)

ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-1.1, 1.1)
ax.set_zlim(-1.1, 1.1)

ax.view_init(elev=30, azim=15)

plt.show()

if True:
    plt.figure()
    plt.plot(degs*RA, degs*dec, '.k')
    plt.xlim(-180, 180)
    plt.ylim(-90, 90)
    plt.show()

> **Comments:**
> **username_0**: Thanks a lot, so you are sure about the dec = np.arcsin(2.*(ran2-0.5)) ? It seems to be good but I don't understand very well why arcsin and not arcos for instance.  
> **username_5**: Very nice! But when I try to rotate it only moves a little bit in the direction of my movement, after which I have to click/drag again for another little bit of movement. Aren't you supposed to be able to drag the plot around smoothly? Is it just me?  
> **username_4**: @pela it works nicely in my IDLE installation which I love, but the 3D plot response in my anaconda installation does not work well at all. This would be a good question for stackoverflow! You'd just edit the script down to the [Minimal, Complete, and Verifiable example](https://stackoverflow.com/help/mcve) before posting. You could link here also, but make sure to get rid of as much as possible.  
> **username_5**: @uhoh: Interesting, I have Anaconda install as well. Good idea with the Stackoverflow, thanks!  

---

### Answer by username_2 on 2017-08-29. Score: 8

Random points on the surface of a sphere can be generated by allowing the azimuthal angle $\phi$ to take a uniformly distributed random value between 0 and $2\pi$. To convert this to RA in degrees you multiply by $180/\pi$. To convert to hours, minutes and seconds you divide the $\phi$ in degrees by 15, which gives the hours, divide the remainder by 60 which gives the minutes and then divide the remainder of that by 60 to give the seconds.

The declination angle $\delta$ is not uniformly distributed, because the area covered at an angle $\delta$ goes as $\cos \delta$.

Note the potential confusion here between declination and the polar angle $\theta$, which is conventionally measured *down* from the pole. The area of a strip of width $\Delta \theta$ at a given $\theta$ on a unit sphere is $$\Delta A = 2\pi \sin \theta\ \Delta \theta = 2\pi \cos (\pi/2 - \delta)\ \Delta \delta = 2\pi \cos \delta\ \Delta \delta$$

To translate this into something we can use, we say that the probability of finding a star in a small range of declinations
$$ dP \propto \cos \delta\ d\delta,$$
so the cumulative probability of finding a value of $\delta$ between $-90^{\circ}$ and $\delta$ is
$$ P(\delta) = \frac{\int^{\delta}\_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}}{\int^{+90^{\circ}}\_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}},$$
where the denominator makes sure the proportionality is correctly normalised. From this, we obtain
$$P(\delta) = \frac{1}{2}(1 + \sin \delta) \tag\*{(1)}$$

The procedure for then obtaining a random declination is to assign a uniform random number to the cumulative probability $P$ between 0 and 1. Then from equation (1) we have
$$ \sin \delta = 2P -1$$
$$\delta = \sin^{-1}(2P-1),$$
where $P$ is a random number between 0 and 1.

> **Comments:**
> **username_2**: @barrycarter Now sorted out - see edit. My original statement was correct.  
> **username_3**: I've just done a Sage / Numpy version that you can rotate, pan & zoom. I'll post the link in the next comment. (Sorry that the code's a bit cryptic, I had to trim it down to make it small enough to fit).  
> **username_3**: [Stars](https://sagecell.sagemath.org/?z=eJxtUstu2zAQvPMrFukhpMvYstwARQABOSXXAOndYCXaJiqRBB9una_vrOSiLtCLsDua2R3tyE0xpEK-TvFCJpOPIvkjdSjWyfghTOvBHkwdyx64dL7IL61SIjpwolt7qcRQTmi2j7QCQhvafm2EgIp-XT5kMpoG26snQZRA60OW3HNrS02e0oqxZJRGmZ1fSi5mnqiazhCeTZL3VZ_vlaiwBRBP4BLvG43NqxYyedb0EN2mZWTTwmhkoyaZyZbk-n0cQ9kN4p3RKNkhz4fy71BNIZrelUvXrLdKvIKa6ySvdP4g5jNT0SEkSoacp4xzHa28WsFHl5Ni8ed_1HW5xrJukaO_0UN24__PGPGMy9tk-jIfdi9zMSl3jw3WfTfj2H1L1Wo6JjfM5XJuM2e0ahEMwrvmuUj5_ry4Q-i49H8p9EBb5r3xmAADu0Hit-jDWCe_B6f_IW8jnlP7sN1OI-UxpO7uU9M0dzzCHWabT5iFe7wvCLu9Iq9A8in8lG-aDpxV92LGbNVvRxXJQQ==&lang=sage)  
> **username_4**: @ProfRob Bravo! Thanks for taking the time towork this out so explicitly. Spherical coordinates tend to start with $\theta = 0$ at the top while latitude and declination put it in the middle, so I always get confused and just start trying $\sin^{-1}$ and $\cos^{-1}$ until something works.  

---

### Answer by username_4 on 2022-03-22. Score: 5

Supplemental answer for future readers here who are more interested in the "uniformness" than the "randomness" or "realism" of the distribution.
Here are two Stack Overflow answers that are similar, but the offset in the sunflower method can be tweaked to optimize the "uniformness" of the nearly uniform distribution.

https://stackoverflow.com/a/44164075/3904031
https://stackoverflow.com/a/26127012/3904031

The main reason that the patterns look different is that the "poles" where the spiral begins and ends is on the side for the "fibonacci_sphere" and the top for the "sunflower_on_sphere".

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

def sunflower_on_sphere(n=1000):
    # https://stackoverflow.com/a/44164075/3904031

    indices = np.arange(n) + 0.5

    phi = np.arccos(1 - 2 * indices / n)
    theta = np.pi * (1 + np.sqrt(5)) * indices

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

    return np.vstack([x, y, z])

def fibonacci_sphere(n=1000):
    # https://stackoverflow.com/a/26127012/3904031
    
    phi = np.pi * (3. - np.sqrt(5.))  # golden angle in radians

    theta = phi * np.arange(n)  # golden angle increment

    y = np.linspace(1, -1, n)  # y goes from 1 to -1
    
    r = np.sqrt(1 - y**2)  # radius at y
    
    x, z = [r*f(theta) for f in (np.cos, np.sin)]

    points = np.vstack([x, y, z])

    return points

n = 2000

xf, yf, zf = fibonacci_sphere(n)
xs, ys, zs = sunflower_on_sphere(n)

fig = plt.figure(figsize=[14, 7.5])

axf  = fig.add_subplot(1, 2, 1, projection='3d', proj_type = 'ortho')
axs  = fig.add_subplot(1, 2, 2, projection='3d', proj_type = 'ortho')

axf.scatter(xf, yf, zf, marker='.', c=yf, cmap='gray')
axs.scatter(xs, ys, zs, marker='.', c=xs, cmap='gray')

for ax in (axf, axs):
    ax.set_xlim(-1.1, 1.1)
    ax.set_ylim(-1.1, 1.1)
    ax.set_zlim(-1.1, 1.1)
    ax.set_box_aspect([1,1,1])
    ax.set_box_aspect([1,1,1]) # https://stackoverflow.com/a/68242226/3904031

axf.view_init(6, -60)
axs.view_init(6, -150)
axf.set_title('fibonacci_sphere(n='+str(n)+')')
axs.set_title('sunflower_on_sphere(n='+str(n)+')')
plt.show()


---

### Answer by username_5 on 2017-08-29. Score: 1

This is really more computer science or mathematics than astronomy, but here's a function random_point_on_unit_sphere()  that creates random points on a unit 3D sphere, which is then used by another function random_point_on_sky() to convert to RA and dec values, using astropy. Finally, the function print_random_star_coords() print a number of RA,dec pairs.
import numpy as np
from astropy.coordinates import SkyCoord
from astropy import units as u

def random_point_on_unit_sphere():
    while True:
        R   = np.random.rand(3) #Random point in box
        R   = 2*R - 1
        rsq = sum(R**2)
        if rsq < 1: break       #Use r only if |r|<1 in order not to favor corners of box
    return R / np.sqrt(rsq)     #Normalize to unit vector

def random_point_on_sky():
    p     = random_point_on_unit_sphere()
    r     = np.linalg.norm(p)
    theta = 90 - (np.arccos(p[2] / r)    / np.pi * 180)            #theta and phi values in degrees
    phi   =       np.arctan(p[1] / p[0]) / np.pi * 180
    c     = SkyCoord(ra=phi, dec=theta, unit=(u.degree, u.degree)) #Create coordinate
    return c.ra.hms, c.dec.deg                                     #Many different formats are possible, e.g c.ra.hour for decimal hour values

def print_random_star_coords(nstars):
    for n in range(nstars):
        RA,dec = random_point_on_sky()
        print(RA,dec)


---

