¿Cómo calcular el ángulo de cono entre dos satélites dados sus ángulos de observación?

nota: Esta respuesta aborda la pregunta directamente:

¿Cómo calcular el ángulo de cono entre dos satélites dados sus ángulos de observación?

Si necesita usar los ángulos de mirada, esta es una buena manera de hacerlo. Esta mejor respuesta le explica al OP que si está usando Skyfield, no debe usar los ángulos de mirada, sino usar las coordenadas en su forma original.

El coseno del ángulo entre dos vectores viene dado por el producto escalar de sus normas.

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

Pero si bajas la distancia y solo usas$az, el$obtienes vectores normalizados automáticamente:

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

¡Creo que esto no es diferente de lo que @AdamTrhon ya ha descrito en esta respuesta ! Puede usar la trigonometría esférica y posiblemente la ley de los cosenos, pero a veces esas expresiones conducen a errores de cálculo debido a cosas como la resta de números casi iguales y la división por casi cero, mientras que de esta manera, trabajando en cartesiano tanto como sea posible, en al menos de la manera que se muestra aquí, no hay posibilidades de que eso suceda.

En Python sería algo como:

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

Entonces, si descarga TLE para tres satélites TDRS más la ISS y calcula el ángulo del cono entre los pares de TDRS, y entre la ISS y cada TDRS, obtendrá algo como lo que se muestra a continuación.

Los objetos también tienen almacenadas sus posiciones geocéntricas, por lo que puede hacer un mapa 3D como se muestra en este ejemplo .

---

EDITAR: lo he usado WhiteSands.at(times).observe(sat.ICRF).apparent().altaz()y esta sería la forma recomendada de obtener la posición óptica aparente. El .apparent()método incluye una variedad de efectos, que incluyen el retraso del tiempo de luz, la aberración astronómica e incluso... espera... efectos gravitatorios de cuerpos masivos que podrían desviar el camino, así como la refracción atmosférica. Puede leer más sobre esto en la documentación en 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()

Aquí está el enfoque manual:

1. Configurar el sistema de coordenadas ortogonales:
   1. La unidad es 1Km (pero no importa mucho).
   2. El origen está en la estación terrestre.
   3. el eje x apunta a 0° de acimut, el eje y a 90°, el eje z hacia arriba.
   4. Hecho: el plano xy es el plano tangente de la Tierra.
2. Calcule los vectores unitarios apuntando de cada satélite. Cuidado con los radianes/grados.
   1. La coordenada x del vector unitario es el coseno del acimut.
   2. La coordenada y del vector unitario es el seno del acimut.
   3. Hecho: El vector ahora apunta al satélite, pero solo en el plano xy.
   4. La coordenada z del vector es tangente de elevación.
   5. Realidad: Ahora el vector apunta al satélite en el espacio 3D, pero ya no es un vector unitario.
   6. Normalícelo: calcule su longitud y divida cada coordenada por esta longitud.
3. Calcule el ángulo del cono:
   1. Calcular el producto escalar de los vectores unitarios
   2. Arccos del producto escalar es el ángulo del cono.

Calcular el ángulo entre dos vectores será difícil si primero transformas sus coordenadas x, y y z en ángulos, porque luego tendrás que sumergirte en las fórmulas de la trigonometría esférica. Skyfield considera de forma nativa que todas las posiciones son vectores x, y, z y, a menudo, es más fácil de calcular si los deja como objetos de "posición" hasta que esté listo para mostrar los resultados.

Si calcula las posiciones de dos satélites en relación con un observador de la Tierra, puede obtener el ángulo entre los satélites utilizando el .separation_from()método de Skyfield que llevan los objetos de posición:

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