计算从卫星到地面的矢量偏离最低点的角度

Question

我想编写 Python 代码来指定具有开普勒元素的卫星轨道，指定地球上的一个点的纬度、经度和高度，指定时间，并计算角度两个向量之间：
- 从卫星到地球上指定点的向量
- 从卫星到地球中心的向量。

我知道我可以使用 poliastro 定义轨道并将其传播到指定时间。困难的部分是在同一坐标系中表示卫星和地球点。

Answer 1

poliastro 目前没有指定坐标系。他们聊天室里有人告诉我地球轨道在 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

计算从卫星到地面的矢量偏离最低点的角度

Compute off-nadir angle of vector from satellite to ground

python

satellite

orbit