Erzeugen Sie eine gleichmäßige Verteilung am Himmel

Durch Zulassen des Azimutwinkels können zufällige Punkte auf der Oberfläche einer Kugel erzeugt werden$\phi$einen gleichmäßig verteilten Zufallswert zwischen 0 und zu nehmen$2\pi$. Um dies in Grad in RA umzuwandeln, multiplizieren Sie mit$180/\pi$. Um in Stunden, Minuten und Sekunden umzurechnen, teilen Sie die$\phi$in Grad durch 15, was die Stunden ergibt, den Rest durch 60 teilen, was die Minuten ergibt, und dann den Rest durch 60 teilen, um die Sekunden zu erhalten.

Der Deklinationswinkel$\delta$nicht gleichmäßig verteilt ist, weil die Fläche schräg abgedeckt wird$\delta$geht wie$\cos \delta$.

Beachten Sie hier die mögliche Verwechslung zwischen Deklination und Polarwinkel$\theta$, die herkömmlich von der Stange nach unten gemessen wird. Die Fläche eines Streifens der Breite$\Delta \theta$bei gegeben$\theta$auf einer Einheitskugel ist

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

Um dies in etwas zu übersetzen, das wir verwenden können, sagen wir, dass die Wahrscheinlichkeit, einen Stern in einem kleinen Bereich von Deklinationen zu finden

$$dP \propto \cos \delta\ d\delta,$$ also die kumulative Wahrscheinlichkeit, einen Wert von zu finden $\delta$zwischen $-90^{\circ}$Und $\delta$Ist $$P(\delta) = \frac{\int^{\delta}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}}{\int^{+90^{\circ}}_{-90^{\circ}}\cos \delta^{\prime}\ d\delta^{\prime}},$$ wobei der Nenner dafür sorgt, dass die Proportionalität korrekt normiert ist. Daraus erhalten wir $$P(\delta) = \frac{1}{2}(1 + \sin \delta) \tag*{(1)}$$

Das Verfahren, um dann eine zufällige Deklination zu erhalten, besteht darin, der kumulativen Wahrscheinlichkeit eine einheitliche Zufallszahl zuzuordnen$P$zwischen 0 und 1. Dann haben wir aus Gleichung (1).

$$\sin \delta = 2P -1$$ $$\delta = \sin^{-1}(2P-1),$$ Wo $P$ist eine Zufallszahl zwischen 0 und 1.

Hier ist mehr Python, als Sie mit einem Teleskop erschüttern können. Ich habe gerade den Algorithmus von @ProfRob verwendet . Dies ist nur ein Python-Skript, die eigentliche Antwort auf die Frage ist die Antwort von @ProfRob und ich habe es gerade geschrieben. Die Mathematik hinter der Erzeugung statistischer Gleichverteilungen wird dort auch sehr schön erklärt.

Python befindet sich unter den Plots. Sie können auf einem XY-Diagramm sehen, dass sie oben und unten dünner werden.

Die 3D-Darstellung ist live, Sie können Ihren Cursor verwenden und die Kugel bewegen. Natürlich nicht HIER :-), sondern in Ihrem eigenen Python-Fenster, wenn Sie das Skript ausführen.

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

Ergänzende Antwort für zukünftige Leser hier, die sich mehr für die "Gleichförmigkeit" als für die "Zufälligkeit" oder den "Realismus" der Verteilung interessieren.

Hier sind zwei Stack Overflow-Antworten, die ähnlich sind, aber der Versatz in der Sonnenblumenmethode kann optimiert werden, um die "Gleichmäßigkeit" der nahezu gleichmäßigen Verteilung zu optimieren.

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

Der Hauptgrund dafür, dass die Muster unterschiedlich aussehen, ist, dass die "Pole", an denen die Spirale beginnt und endet, bei der "fibonacci_sphere" auf der Seite und bei der "sunflower_on_sphere" oben liegen.

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

Das ist wirklich mehr Informatik oder Mathematik als Astronomie, aber hier ist eine Funktion random_point_on_unit_sphere(), die zufällige Punkte auf einer 3D-Einheitskugel erstellt, die dann von einer anderen Funktion verwendet wird, random_point_on_sky()um sie in RA- und Dez-Werte umzuwandeln, wobei Astropie verwendet wird . Schließlich gibt die Funktion print_random_star_coords()eine Reihe von RA,dec-Paaren aus.

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)