Вычислить угол отклонения вектора от спутника к земле

Question

Я хотел бы написать код Python, чтобы указать орбиту спутника с кеплеровскими элементами, указать точку на Земле с широтой, долготой и высотой, указать время и вычислить угол между двумя векторами:
- вектор от спутника до указанной точки на Земле
- вектор от спутника до центра Земли.

Я знаю, что могу использовать полиастро для определения орбиты и распространения ее на указанное время. Трудная часть представляет спутник и точку Земли в одной и той же системе координат.

Answer 1

В настоящее время полиастро не определяет систему координат. Кто-то в их чате сказал мне, что земные орбиты находятся в GCRS. astropy может преобразовать GCRS в ITRS, которая представляет собой ориентированную на Землю систему координат:

import math
import numpy as np
from astropy import units as u
from astropy.time import Time
from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from astropy.coordinates import SkyCoord

def lla2ecef(lat, lon, alt):
    """ Convert lat/lon/alt to cartesian position in ECEF frame.

    Origin is center of Earth. +x axis goes through lat/lon (0, 0).
    +y axis goes through lat/lon (0, 90). +z axis goes through North Pole.

    lat: number, geodetic latitude in degrees
    lon: number, longitude in degrees
    alt: number, altitude above WGS84 ellipsoid, in km

    Returns: tuple (x, y, z) coordinates, in km.

    Source: "Department of Defense World Geodetic System 1984"
    Page 4-4
    National Imagery and Mapping Agency
    Last updated June, 2004
    NIMA TR8350.2
    """

    lon = lon * math.pi/180.0  # Convert to radians
    lat = lat * math.pi/180.0  # Convert to radians

    # WGS84 ellipsoid constants:
    a = 6378.137 #equatorial radius, in km
    e = 8.1819190842622e-2

    # intermediate calculation: prime vertical radius of curvature
    N =  a/math.sqrt(1 - e**2*math.sin(lat)**2)

    #results
    x = (N + alt)*math.cos(lat)*math.cos(lon)
    y = (N + alt)*math.cos(lat)*math.sin(lon)
    z = ((1 - e**2)*N + alt)*math.sin(lat)

    return (x, y, z)

epoch = Time(2018, format='decimalyear', scale='tai') #01-01-2018 00:00:00, in TAI
propagation_time = 9000 #seconds
semi_major_axis = 10000 #km
eccentricity = 0.1
inclination = 50 #deg
raan = 70 #deg
arg_perigee = 60 #deg
true_anomaly = 80 #deg
orbit = Orbit.from_classical(
    Earth, 
    semi_major_axis*u.km, 
    eccentricity*u.one,  
    inclination*u.deg, raan*u.deg, 
    arg_perigee*u.deg, 
    true_anomaly*u.deg, 
    epoch)
propagated_orbit = orbit.propagate(propagation_time*u.s)
pos_gcrs = propagated_orbit.state.r
sky_gcrs = SkyCoord(
    representation_type='cartesian', 
    x=pos_gcrs[0], y=pos_gcrs[1], z=pos_gcrs[2],
    frame='gcrs',
    obstime=(epoch + propagation_time*u.s))
pos_ecef = sky_gcrs.transform_to('itrs')
pos_s = np.array((pos_ecef.x.to(u.km).value, 
                  pos_ecef.y.to(u.km).value, 
                  pos_ecef.z.to(u.km).value))

lat = 40 #deg
lon = 50 #deg
alt = 0.06 #km
pos_t = np.array(lla2ecef(lat, lon, alt))

#Compute angle at satellite between target and center of earth
v1 = pos_t - pos_s
v2 = -pos_s
angle = math.acos(np.dot(v1, v2)/(np.linalg.norm(v1)*np.linalg.norm(v2)))

#convert to degrees
angle = angle*180.0/math.pi