10
$\begingroup$

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!

$\endgroup$
4
  • $\begingroup$ What about DEC=Arccos(np.random.uniform(-90,90))? $\endgroup$
    – J. Chomel
    Aug 29, 2017 at 9:08
  • $\begingroup$ 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. $\endgroup$
    – Richard
    Aug 29, 2017 at 9:44
  • $\begingroup$ 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 ? $\endgroup$
    – Richard
    Aug 29, 2017 at 10:10
  • $\begingroup$ 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. $\endgroup$
    – user21
    Aug 29, 2017 at 12:03

4 Answers 4

8
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – user21
    Aug 29, 2017 at 11:59
  • $\begingroup$ @barrycarter Now sorted out - see edit. My original statement was correct. $\endgroup$
    – ProfRob
    Aug 29, 2017 at 14:57
  • $\begingroup$ 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). $\endgroup$
    – PM 2Ring
    Jul 10, 2021 at 16:25
  • $\begingroup$ Stars $\endgroup$
    – PM 2Ring
    Jul 10, 2021 at 16:25
  • $\begingroup$ @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. $\endgroup$
    – uhoh
    Apr 6, 2022 at 23:39
6
$\begingroup$

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.

enter image description here

enter image description here

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()
$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – Richard
    Aug 29, 2017 at 13:23
  • $\begingroup$ 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? $\endgroup$
    – pela
    Aug 29, 2017 at 13:53
  • 1
    $\begingroup$ OK, I get it: when you integrate the cos, you get -sin, so arcsin is apparently the right one to use. $\endgroup$
    – user21
    Aug 29, 2017 at 15:25
  • 1
    $\begingroup$ @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. $\endgroup$
    – uhoh
    Aug 29, 2017 at 15:35
  • 1
    $\begingroup$ @uhoh: Interesting, I have Anaconda install as well. Good idea with the Stackoverflow, thanks! $\endgroup$
    – pela
    Aug 29, 2017 at 18:06
5
$\begingroup$

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

two nearly uniform distributions of points on a 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()
$\endgroup$
1
$\begingroup$

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)
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .