Generate a uniform distribution on the sky

Question

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!

Answer 1

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.

Answer 2

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

Answer 3

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 4

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)