两个Topos（）位置之间的路径：确定到达给定高度的纬度和经度

Question

我是轨道力学领域的新手，目前正在努力解决以下问题，这应该很容易用Skyfield解决，但我有点被所有不同的坐标系和它们之间的转换所淹没

我在地球上有一个{}位置和一个低轨卫星的{}位置。我正在考虑他们之间的视线。我想确定沿着这条路径的纬度和经度，它与大气的特定层相交的位置

一个例子是中层和一个基于纬度和经度给出的关于其在大约100公里处的属性的现有数据集。交叉点将使我更好地理解这些属性在与卫星通信中的相互作用

我试着直接用Skyfield来做，但最终得到的是一个Apparent对象，我无法转换回地球上的纬度和经度。首先，我用三角法确定了从地球到达到100公里高度的点的距离

然后，我选择了地球上的位置，使用不变的仰角、方位角来保持路径的方向，最后加上计算出的距离来到达这个位置。我想我需要得到一个Geocentric对象来使用subpoint()，以便获得这个位置所需的纬度和经度

这就是我到目前为止所做的：

from skyfield.api import load, Distance
from skyfield.toposlib import Topos
import numpy as np

ts = load.timescale()

earth_position = Topos('52.230039 N', '4.842402 E', elevation_m=10)
space_position = Topos('51.526200 N', '5.347795 E', elevation_m=625 * 1000)

difference = (space_position - earth_position).at(ts.now()).altaz()

distance_to_height = 100 / np.sin(difference[0].radians)

position = earth_position.at(ts.now()).from_altaz(alt_degrees=difference[0].degrees, az_degrees=difference[1].degrees, distance=Distance(km=distance_to_height))

我已经阅读了多次文档，偶然发现了frame_latlon(frame)通用ICRF对象，但不确定如何进一步进行

尝试将纬度和经度完全三角化也不会产生预期的结果

不幸的是，我没有任何验证结果可以用来更容易地解决这个问题。从三角角度再想象一下，很明显，卫星位置的高度增加会使交叉点的纬度、经度更接近地球上的位置。降低高度将使该交叉点靠近卫星

Answer 1

这是一个有趣的问题，Skyfield的API并没有提供简单的方法来询问；如果您可以通过了解视线与特定高度的交点来概述将要解决的更大问题，那么有可能为解决同一问题的未来用户编写解决该问题的例行程序

同时：

1. 为了让脚本运行，我必须从api导入Distance
2. 名称dis无法识别，因此我将其替换为distance_to_height，希望它是预期的名称
3. 调用ts.now()是指每次调用的日期和时间略有不同。虽然脚本运行速度如此之快，可能并不重要，但为了清楚起见，我只在脚本开始时调用now()一次，这也比反复调用要快一些。（实际上，在这种情况下，速度要快得多，因为旋转矩阵只计算一次，而不必为每个单独的时间对象重复计算，但这是一个不容易看到的隐藏细节。）
4. 我怀疑你的几何学有问题：我想只有在地球是平的情况下100 / sin()机动才有效？但也许你总是在处理接近头顶的卫星，所以误差是可控的？（或者可能是我错误地想象了几何体；如果数学事实上是正确的，请随意提供图表。）
5. 为了可读性，我将给出altaz()的组件名称，而不是数字

有了这些调整， 我认为答案是您需要手动构造一个Geocentric 通过将观察者的位置相加来确定位置 以及您创建的相对向量 在观察者和视线上100公里点之间。 必须采取像这样的手动步骤意味着一个可能的领域 Skyfield可以改进的地方。 下面是它在代码中的外观：

from skyfield.api import load, Distance
from skyfield.positionlib import Geocentric
from skyfield.toposlib import Topos
import numpy as np

ts = load.timescale()
t = ts.now()

earth_position = Topos('52.230039 N', '4.842402 E', elevation_m=10)
space_position = Topos('51.526200 N', '5.347795 E', elevation_m=625 * 1000)

alt, az, distance = (space_position - earth_position).at(t).altaz()

distance_to_height = 100 / np.sin(alt.radians)

e = earth_position.at(t)
p = e.from_altaz(alt_degrees=alt.degrees, az_degrees=az.degrees, distance=Distance(km=distance_to_height))

g = Geocentric(e.position.au + p.position.au, t=t)
s = g.subpoint()
print(s)
print(s.elevation.km, '<- warning: 100/sin() did not produce exactly 100')

我看到的结果是：

Topos 52deg 06' 30.0" N 04deg 55' 51.7" E
100.02752954478532 <- warning: 100/sin() did not produce exactly 100

为了未来,， 我在Skyfield TODO.rst文件中添加了一些想法 这两者结合起来可能会朝着一种更惯用的方式发展 在将来执行这种计算- 尽管我怀疑，除此之外，还有一些步骤是必要的：

https://github.com/skyfielders/python-skyfield/commit/ba1172a0ccfef84473436d9d7b8a7d7011344cbd