从卫星到地面的矢量离地角的计算

Question

我想编写Python代码来指定一个卫星的轨道和Keplerian，指定地球上一个纬度、经度和高度的点，指定一个时间，并计算两个矢量之间的角度：

• 从卫星到地球上指定点的矢量
• 从卫星到地球中心的矢量。

我知道我可以用小儿麻痹症来定义轨道并将其传播到指定的时间。硬部分是表示卫星和地球点在同一个坐标系中。

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