Generar una distribución uniforme en el cielo

Se pueden generar puntos aleatorios en la superficie de una esfera permitiendo que el ángulo azimutal$\phi$para tomar un valor aleatorio uniformemente distribuido entre 0 y$2\pi$. Para convertir esto a RA en grados, multiplique por$180/\pi$. Para convertir a horas, minutos y segundos se divide el$\phi$en grados por 15, lo que da las horas, divide el resto por 60 que da los minutos y luego divide el resto por 60 para dar los segundos.

El ángulo de declinación$\delta$no se distribuye uniformemente, porque el área cubierta en un ángulo$\delta$va como$\cos \delta$.

Tenga en cuenta la posible confusión aquí entre la declinación y el ángulo polar$\theta$, que se mide convencionalmente hacia abajo desde el poste. El área de una franja de ancho$\Delta \theta$en un dado$\theta$en una esfera unitaria es

$$\Delta A = 2\pi \sin \theta\ \Delta \theta = 2\pi \cos (\pi/2 - \delta)\ \Delta \delta = 2\pi \cos \delta\ \Delta \delta$$

Para traducir esto en algo que podamos usar, decimos que la probabilidad de encontrar una estrella en un pequeño rango de declinaciones

$$dP \propto \cos \delta\ d\delta,$$ entonces la probabilidad acumulada de encontrar un valor de $\delta$entre $-90^{\circ}$y $\delta$es $$P(\delta) = \frac{\int^{\delta}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}}{\int^{+90^{\circ}}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}},$$ donde el denominador asegura que la proporcionalidad esté correctamente normalizada. De esto, obtenemos $$P(\delta) = \frac{1}{2}(1 + \sin \delta) \tag*{(1)}$$

El procedimiento para obtener entonces una declinación aleatoria es asignar un número aleatorio uniforme a la probabilidad acumulada$P$entre 0 y 1. Entonces de la ecuación (1) tenemos

$$\sin \delta = 2P -1$$ $$\delta = \sin^{-1}(2P-1),$$ dónde $P$es un número aleatorio entre 0 y 1.

Aquí hay más pitones de los que puedes sacudir con un telescopio. Acabo de usar el algoritmo de @ProfRob . Esta es solo una secuencia de comandos de Python, la respuesta real a la pregunta es la respuesta de @ProfRob y la acabo de escribir. Las matemáticas detrás de la generación de distribuciones estadísticamente uniformes también se explican muy bien allí.

Python está debajo de las parcelas. Puede ver en un gráfico XY que se adelgazan en la parte superior e inferior.

La trama 3D está en vivo, puede usar el cursor y mover la esfera. Por supuesto, no AQUÍ :-), sino en su propia ventana de Python si ejecuta el 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()

Respuesta complementaria para futuros lectores aquí que estén más interesados ​​​​en la "uniformidad" que en la "aleatoriedad" o el "realismo" de la distribución.

Aquí hay dos respuestas de desbordamiento de pila que son similares, pero el desplazamiento en el método del girasol se puede modificar para optimizar la "uniformidad" de la distribución casi uniforme.

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

La razón principal por la que los patrones se ven diferentes es que los "polos" donde comienza y termina la espiral están en el costado de la "esfera_fibonacci" y en la parte superior para el "girasol_en_esfera".

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

Esto es más ciencia informática o matemáticas que astronomía, pero aquí hay una función random_point_on_unit_sphere()que crea puntos aleatorios en una esfera 3D unitaria, que luego es utilizada por otra función random_point_on_sky()para convertir a valores RA y dec, usando astropy . Finalmente, la función print_random_star_coords()imprime un número de pares RA,dec.

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)