Как рассчитать угол конуса между двумя спутниками, учитывая их углы обзора?

примечание: этот ответ напрямую касается вопроса:

Как рассчитать угол конуса между двумя спутниками, учитывая их углы обзора?

Если вам нужно использовать углы обзора, это хороший способ сделать это. Этот лучший ответ объясняет ОП, что если вы используете Skyfield, вам не следует использовать углы обзора, а вместо этого использовать координаты в их исходной форме.

Косинус угла между двумя векторами определяется скалярным произведением их норм.

$$\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)$$

Но если вы сократите расстояние и просто используете$az, el$вы автоматически получаете нормализованные векторы:

$$\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}}$$

Я думаю, что это ничем не отличается от того, что @AdamTrhon уже описал в этом ответе ! Вы можете использовать сферическую тригонометрию и, возможно, закон косинусов, но иногда эти выражения приводят к вычислительным ошибкам из-за таких вещей, как вычитание почти равных чисел и деление почти на ноль, тогда как таким образом - работая в декартовой системе, насколько это возможно, в способ, показанный здесь, по крайней мере, - нет никаких шансов, что это произойдет.

В Python это будет примерно так:

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

Итак, если вы загрузите TLE для трех спутников TDRS плюс ISS и вычислите угол конуса между парами TDRS, а также между ISS и каждым TDRS, вы получите что-то вроде того, что показано ниже.

Объекты также сохраняют свое геоцентрическое положение, поэтому вы можете создать 3D-карту, как показано в этом примере .

---

РЕДАКТИРОВАТЬ: я использовал, WhiteSands.at(times).observe(sat.ICRF).apparent().altaz()и это был бы рекомендуемый способ получить видимое оптическое положение. Метод .apparent()включает множество эффектов, в том числе временную задержку света, астрономическую аберрацию и даже... подождите... гравитационные эффекты массивных тел, которые могут отклонить траекторию, а также атмосферную рефракцию. Подробнее об этом можно прочитать в документации по адресу http://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.Astrometric.apparent .

---

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

Вот ручной подход:

1. Настройте ортогональную систему координат:
   1. Единица 1 км (но это не имеет большого значения).
   2. Происхождение находится на наземной станции.
   3. ось x указывает на азимут 0°, ось y на 90°, ось z строго вверх.
   4. Факт: плоскость xy является касательной к Земле.
2. Вычислите единичные векторы, указывающие на каждый спутник. Остерегайтесь радианов/градусов.
   1. Координата x единичного вектора представляет собой косинус азимута.
   2. Координата y единичного вектора представляет собой синус азимута.
   3. Факт: вектор теперь указывает на спутник, но только в плоскости xy.
   4. Z-координата вектора является касательной к высоте.
   5. Факт: теперь вектор указывает на спутник в трехмерном пространстве, но это уже не единичный вектор.
   6. Нормализуйте его: вычислите его длину и разделите каждую координату на эту длину.
3. Вычислите угол конуса:
   1. Вычислить скалярное произведение единичных векторов
   2. Arccos скалярного произведения представляет собой угол конуса.

Вычисление угла между двумя векторами будет затруднено, если вы сначала преобразуете их координаты x, y и z в углы, потому что тогда вам нужно погрузиться в формулы сферической тригонометрии. Skyfield изначально рассматривает все позиции как векторы x, y, z, и часто их проще вычислить, если вы оставите их как объекты «позиции», пока не будете готовы отобразить результаты.

Если вы вычисляете положение двух спутников относительно наблюдателя Земли, вы можете получить угол между спутниками, используя .separation_from()метод Скайфилда, который несут объекты положения:

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