# Generate an uniform distribution on the sky

How to generate an 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))? – J. Chomel Aug 29 '17 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. – Richard Aug 29 '17 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 ? – Richard Aug 29 '17 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 '17 at 12:03

Heres more python than you can shake a telescope at. I just used @RobJeffries' algorithm. This is just a python script, the real answer to the question is @RobJeffries' 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. – Richard Aug 29 '17 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 '17 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 '17 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 '17 at 15:35
• @Richard I'd strongly recommend accepting RobJeffries' answer instead of this one. This is just a script, it does not explain very much. The other answer really answers your question about how to do it in a way that can be generalized. – uhoh Aug 29 '17 at 16:40

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 '17 at 11:59
• @barrycarter there's a line with your name on it in my script. $dec = sin^{-1}(random)$ is correct; Rob Jeffries is not wrong. – uhoh Aug 29 '17 at 12:04
• @barrycarter Now sorted out - see edit. My original statement was correct. – Rob Jeffries Aug 29 '17 at 14:57
• @RobJeffries Bravo! Thanks for taking the time to work 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 Aug 29 '17 at 15:32

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 / r)    / np.pi * 180)            #theta and phi values in degrees
phi   =       np.arctan(p / p) / 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)