# Generate a uniform distribution on the sky

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)...

• What about DEC=Arccos(np.random.uniform(-90,90))? Aug 29, 2017 at 9:08
• 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. Aug 29, 2017 at 9:44
• 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 ? Aug 29, 2017 at 10:10
• I'd go with arccos because it makes more sense (length of a declination line is proportional to cos, not sin), but I'm guessing arcsin yields the same distribution.
– user21
Aug 29, 2017 at 12:03

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.

• I think you mean the area covered at angle delta is cos(delta), not sin(delta), since the area (length) covered by the celestial equator at 0 degrees is the highest value, not the lowest value.
– user21
Aug 29, 2017 at 11:59
• @barrycarter Now sorted out - see edit. My original statement was correct. Aug 29, 2017 at 14:57
• 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). Jul 10, 2021 at 16:25
• Stars Jul 10, 2021 at 16:25
• @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
Apr 6 at 23:39

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]]

# 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()

• 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. Aug 29, 2017 at 13:23
• 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
Aug 29, 2017 at 13:53
• OK, I get it: when you integrate the cos, you get -sin, so arcsin is apparently the right one to use.
– user21
Aug 29, 2017 at 15:25
• @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 before posting. You could link here also, but make sure to get rid of as much as possible.
– uhoh
Aug 29, 2017 at 15:35
• @uhoh: Interesting, I have Anaconda install as well. Good idea with the Stackoverflow, thanks!
– pela
Aug 29, 2017 at 18:06

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.

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()


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)