Как создать равномерное распределение звезд на небе? Я хочу создать симуляцию случайных точек, следующих за однородным распределением на части неба (сделав предположение, что у нас одинаковое количество звезд в любом направлении). Проблема в том, как мне это сделать из-за cos(dec) в RA?
То, что я сделал в питоне, это
import numpy as np
RA = np.random.uniform(-180,180)
Dec = np.random.uniform(-90,90)
Но это совершенно неверно, так как у нас будет больше точек на полюсах (если я сделаю много точек, это будет явно неравномерно)...
Заранее большое спасибо!
Случайные точки на поверхности сферы можно генерировать, допуская азимутальный угол чтобы взять равномерно распределенное случайное значение между 0 и . Чтобы преобразовать это в RA в градусах, вы умножаете на . Чтобы преобразовать в часы, минуты и секунды, вы делите в градусах на 15, что дает часы, разделите остаток на 60, чтобы получить минуты, а затем разделите остаток на 60, чтобы получить секунды.
Угол склонения распределяется неравномерно, потому что площадь покрытия под углом идет как .
Обратите внимание на возможную путаницу между склонением и полярным углом. , который обычно измеряется вниз от полюса. Площадь полосы шириной в данный на единичной сфере
Чтобы перевести это в нечто, что мы можем использовать, мы говорим, что вероятность найти звезду в небольшом диапазоне склонений
Процедура получения случайного склонения состоит в том, чтобы присвоить единообразное случайное число кумулятивной вероятности между 0 и 1. Тогда из уравнения (1) имеем
cos(delta)
, а не sin(delta)
, поскольку площадь (длина), покрытая небесным экватором при 0 градусах, является самым высоким значением, а не самым низким значением.Здесь больше питона, чем вы можете потрясти телескопом. Я только что использовал алгоритм @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, которые похожи, но смещение в методе подсолнуха можно настроить, чтобы оптимизировать «равномерность» почти равномерного распределения.
Основная причина того, что паттерны выглядят по-разному, заключается в том, что «полюса», где спираль начинается и заканчивается, находятся сбоку для «сферы Фибоначчи» и сверху для «подсолнуха на_сфере».
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)
Дж. Хомель
DEC=Arccos(np.random.uniform(-90,90))
?Ричард
Ричард
пользователь21
arccos
потому что это имеет больше смысла (длина линии склонения пропорциональнаcos
, а неsin
), но я предполагаю, чтоarcsin
дает такое же распределение.