Создайте равномерное распределение на небе

Случайные точки на поверхности сферы можно генерировать, допуская азимутальный угол$\phi$чтобы взять равномерно распределенное случайное значение между 0 и$2\pi$. Чтобы преобразовать это в RA в градусах, вы умножаете на$180/\pi$. Чтобы преобразовать в часы, минуты и секунды, вы делите$\phi$в градусах на 15, что дает часы, разделите остаток на 60, чтобы получить минуты, а затем разделите остаток на 60, чтобы получить секунды.

Угол склонения$\delta$распределяется неравномерно, потому что площадь покрытия под углом$\delta$идет как$\cos \delta$.

Обратите внимание на возможную путаницу между склонением и полярным углом.$\theta$, который обычно измеряется вниз от полюса. Площадь полосы шириной$\Delta \theta$в данный$\theta$на единичной сфере

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

Чтобы перевести это в нечто, что мы можем использовать, мы говорим, что вероятность найти звезду в небольшом диапазоне склонений

$$dP \propto \cos \delta\ d\delta,$$ поэтому кумулятивная вероятность найти значение $\delta$между $-90^{\circ}$и $\delta$является $$P(\delta) = \frac{\int^{\delta}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}}{\int^{+90^{\circ}}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}},$$ где знаменатель обеспечивает правильную нормализацию пропорциональности. Отсюда получаем $$P(\delta) = \frac{1}{2}(1 + \sin \delta) \tag*{(1)}$$

Процедура получения случайного склонения состоит в том, чтобы присвоить единообразное случайное число кумулятивной вероятности$P$между 0 и 1. Тогда из уравнения (1) имеем

$$\sin \delta = 2P -1$$ $$\delta = \sin^{-1}(2P-1),$$ где $P$случайное число от 0 до 1.

Здесь больше питона, чем вы можете потрясти телескопом. Я только что использовал алгоритм @ProfRob . Это всего лишь скрипт на Python, настоящий ответ на вопрос — это ответ @ProfRob , и я только что написал его. Там также очень хорошо объясняется математика создания статистически однородных распределений.

Python находится ниже графиков. Вы можете видеть на графике XY, что они утончаются вверху и внизу.

3D-сюжет живой, вы можете использовать свой курсор и перемещать сферу. Конечно, не ЗДЕСЬ :-), а в вашем собственном окне Python, если вы запустите скрипт.

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

Дополнительный ответ для будущих читателей, которых больше интересует «однородность», чем «случайность» или «реализм» распределения.

Вот два ответа Stack Overflow, которые похожи, но смещение в методе подсолнуха можно настроить, чтобы оптимизировать «равномерность» почти равномерного распределения.

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

Основная причина того, что паттерны выглядят по-разному, заключается в том, что «полюса», где спираль начинается и заканчивается, находятся сбоку для «сферы Фибоначчи» и сверху для «подсолнуха на_сфере».

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

Это действительно больше информатика или математика, чем астрономия, но вот функция random_point_on_unit_sphere(), которая создает случайные точки на единичной трехмерной сфере, которые затем используются другой функцией random_point_on_sky()для преобразования в значения RA и dec, используя astropy . Наконец, функция print_random_star_coords()выводит количество пар 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)