Вычислить угол отклонения вектора от спутника к земле
Я хотел бы написать код Python, чтобы указать орбиту спутника с кеплеровскими элементами, указать точку на Земле с широтой, долготой и высотой, указать время и вычислить угол между двумя векторами:
- вектор от спутника до указанной точки на Земле
- вектор от спутника до центра Земли.
Я знаю, что могу использовать полиастро для определения орбиты и распространения ее на указанное время. Трудная часть представляет спутник и точку Земли в одной и той же системе координат.
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