Wie berechnet man den Kegelwinkel zwischen zwei Satelliten aufgrund ihrer Blickwinkel?

Hinweis: Diese Antwort adressiert direkt die Frage:

Wie berechnet man den Kegelwinkel zwischen zwei Satelliten aufgrund ihrer Blickwinkel?

Wenn Sie die Blickwinkel verwenden müssen, ist dies eine gute Möglichkeit, dies zu tun. Diese bessere Antwort erklärt dem OP, dass Sie bei Verwendung von Skyfield nicht die Blickwinkel , sondern die Koordinaten in ihrer ursprünglichen Form verwenden sollten.

Der Kosinus des Winkels zwischen zwei Vektoren ist durch das Skalarprodukt ihrer Normen gegeben.

$$\cos(\theta_{12}) = \mathbf{\hat{v}_1} \cdot \mathbf{\hat{v}_2} = \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{|\mathbf{v_1}| \ |\mathbf{v_2}|} = \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{v_1 \ v_2}$$

$$\theta_{12} = \cos^{-1}\left( \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{v_1 \ v_2} \right)$$

Aber wenn Sie die Entfernung fallen lassen und einfach verwenden$az, el$Sie erhalten automatisch einen normalisierten Vektor:

$$\mathbf{\hat{v}_i} = \cos(el_i)\left(\cos(az_i)\mathbf{\hat{x}}+\sin(az_i)\mathbf{\hat{y}}\right) + \sin(el_i)\mathbf{\hat{z}}$$

Ich denke, das ist nicht anders als das, was @AdamTrhon bereits in dieser Antwort beschrieben hat ! Sie können sphärischen Trig und möglicherweise den Kosinussatz verwenden, aber manchmal führen diese Ausdrücke zu Rechenfehlern aufgrund von Dingen wie der Subtraktion von nahezu gleichen Zahlen und der Division durch nahezu Null, während auf diese Weise - so viel wie möglich kartesisch gearbeitet wird, in zumindest der hier gezeigte Weg, - dafür gibt es keine Chance.

In Python wäre das etwa so:

def nvec(elaz):
    (cel, sel), (caz, saz) = [[f(q) for f in (np.cos, np.sin)] for q in elaz]   # parentheses for Py3
    return np.array([cel*caz, cel*saz, sel])

def angle(elaz1, elaz2):
    v1, v2 = [nvec(elaz) for elaz in (elaz1, elaz2)]  # parentheses for Py3
    return np.arccos((v1*v2).sum(axis=0))

Wenn Sie also TLEs für drei TDRS-Satelliten plus die ISS herunterladen und den Kegelwinkel zwischen TDRS-Paaren und zwischen der ISS und jedem TDRS berechnen, erhalten Sie so etwas wie das unten gezeigte.

Die geozentrischen Positionen der Objekte werden ebenfalls gespeichert, sodass Sie eine 3D-Karte erstellen können, wie in diesem Beispiel gezeigt .

---

EDIT: Ich habe verwendet WhiteSands.at(times).observe(sat.ICRF).apparent().altaz()und dies wäre der empfohlene Weg, um die scheinbare optische Position zu erhalten. Die .apparent()Methode umfasst eine Vielzahl von Effekten, darunter die Lichtzeitverzögerung, astronomische Aberration und sogar ... warten Sie darauf ... Gravitationseffekte massiver Körper, die den Weg ablenken könnten, sowie atmosphärische Brechung. Sie können mehr darüber in der Dokumentation unter http://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.Astrometric.apparent nachlesen

---

TLEs = """TDRS 5                  
1 21639U 91054B   18086.36437858  .00000071  00000-0  00000-0 0  9995
2 21639  14.5306  18.4626 0026343 345.4651 144.6288  1.00281508 97593
TDRS 10                 
1 27566U 02055A   18087.12756861  .00000056  00000-0  00000+0 0  9998
2 27566   5.5204  57.1630 0011308 250.2541 109.7547  1.00278469 56111
TDRS 11                 
1 39070U 13004A   18086.87347718  .00000063  00000-0  00000-0 0  9994
2 39070   5.0128 328.6219 0008993 321.0893  38.7221  1.00272889 16583
ISS (ZARYA)             
1 25544U 98067A   18088.22902370  .00003740  00000-0  63642-4 0  9999
2 25544  51.6415  57.3234 0001506 271.5382 195.6957 15.54152785106088"""

lines           = TLEs.splitlines()
names, L1s, L2s = [[x.strip() for x in lines[i::3]] for i in range(3)]
triplets        = zip(names, L1s, L2s)

class Sat(object):
    def __init__(self, name):
        self.name = name

def nvec(elaz):
    (cel, sel), (caz, saz) = [[f(q) for f in (np.cos, np.sin)] for q in elaz]   # parentheses for Py3
    return np.array([cel*caz, cel*saz, sel])

def angle(elaz1, elaz2):
    v1, v2 = [nvec(elaz) for elaz in (elaz1, elaz2)]  # parentheses for Py3
    return np.arccos((v1*v2).sum(axis=0))

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos, EarthSatellite
import itertools

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]  # parentheses for Py3
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
ts   = load.timescale()

data    = load('de421.bsp')
earth   = data['earth']
ts      = load.timescale()

WhiteSands  = earth + Topos(latitude_degrees   =   32.4,
                             longitude_degrees = -106.5,
                             elevation_m       = 1300.0  )

minutes = np.arange(0, 1441, 1)
times   = ts.utc(2018, 3, 29, 0, minutes)

sats = []
for name, L1, L2 in triplets:
    sat = Sat(name)
    sats.append(sat)
    sat.Geo  = EarthSatellite(L1, L2)
    sat.ICRF = earth + EarthSatellite(L1, L2)
    sat.obs  = WhiteSands.at(times).observe(sat.ICRF)
    sat.elaz = [x.radians for x in sat.obs.apparent().altaz()[:2]]
    sat.below = sat.elaz[0] <= 0.
    sat.pos   = sat.Geo.at(times).position.km

ISS   = [sat for sat in sats if 'ISS' in sat.name][0]
TDRSs = [sat for sat in sats if 'TDRS' in sat.name]
TDRSpairs = list(itertools.combinations(TDRSs, 2))

intra_TDRS_cones = []
for pair in TDRSpairs:
    name = ''.join([x.name + ' ' for x in pair])[:-1]
    elaz1, elaz2 = [s.elaz for s in pair]
    cone = angle(elaz1, elaz2)
    intra_TDRS_cones.append((name, cone))

ISS_TDRS_cones = []
for TDRS in TDRSs:
    name = ISS.name + ' ' + TDRS.name
    cone = angle(ISS.elaz, TDRS.elaz)
    ISS_TDRS_cones.append((name, cone))

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    for name, cone in intra_TDRS_cones:
        plt.plot(minutes/60., degs*cone)
    plt.title("intra-TDRS cone angles", fontsize=16)
    plt.xlim(0, 24)
    plt.subplot(2, 1, 2)
    for name, cone in ISS_TDRS_cones:
        plt.plot(minutes/60., degs*cone)
    plt.title("ISS-TDRS cone angles", fontsize=16)
    plt.xlim(0, 24)
    plt.show()

Hier ist der manuelle Ansatz:

1. Orthogonales Koordinatensystem einrichten:
   1. Die Einheit ist 1 km (aber es spielt keine große Rolle).
   2. Ursprung ist in der Bodenstation.
   3. x-Achse zeigt auf 0° Azimut, y-Achse auf 90°, z-Achse senkrecht nach oben.
   4. Tatsache: Die xy-Ebene ist die Tangentialebene der Erde.
2. Berechnen Sie Einheitsvektoren, die auf jeden Satelliten zeigen. Vorsicht bei Bogenmaß/Grad.
   1. Die x-Koordinate des Einheitsvektors ist der Kosinus des Azimuts.
   2. Die y-Koordinate des Einheitsvektors ist der Sinus des Azimuts.
   3. Tatsache: Der Vektor zeigt jetzt auf den Satelliten, aber nur in der xy-Ebene.
   4. Die z-Koordinate des Vektors ist die Tangente der Höhe.
   5. Tatsache: Jetzt zeigt der Vektor auf den Satelliten im 3D-Raum, aber es ist kein Einheitsvektor mehr.
   6. Normalisieren Sie es: Berechnen Sie seine Länge und teilen Sie jede Koordinate durch diese Länge.
3. Kegelwinkel berechnen:
   1. Berechnen Sie das Skalarprodukt der Einheitsvektoren
   2. Arccos des Punktprodukts ist der Kegelwinkel.

Das Berechnen des Winkels zwischen zwei Vektoren wird schwierig, wenn Sie zuerst ihre x-, y- und z-Koordinaten in Winkel umwandeln, weil Sie dann in die Formeln der sphärischen Trigonometrie eintauchen müssen. Skyfield betrachtet alle Positionen nativ als x-, y-, z-Vektoren, und es ist oft einfacher zu berechnen, wenn Sie sie als „Positions“-Objekte belassen, bis Sie bereit sind, Ergebnisse anzuzeigen.

Wenn Sie zwei Satellitenpositionen relativ zu einem Erdbeobachter berechnen, können Sie den Winkel zwischen den Satelliten mithilfe der Skyfield- .separation_from()Methode erhalten, die Positionsobjekte tragen:

from skyfield import api

# Time.
ts = api.load.timescale()
t = ts.utc(2018, 3, 30, 23, 8)

# Satellites.
sats = api.load.tle('https://celestrak.com/NORAD/elements/stations.txt')
s1 = sats['ISS']
s2 = sats['ASTERIA']

# Observe the satellites from a position on the Earth's surface.
usno = api.Topos('38.9215 N', '77.0669 W', elevation_m=92.0)
pos1 = (s1 - usno).at(t)
pos2 = (s2 - usno).at(t)

# How far apart are the satellites in the sky?
print(pos1.separation_from(pos2))