Path between two Topos() locations: determine latitude and longitude where a given altitude is reached

Question

I am quite new to field of orbital mechanics and currently struggling a bit with the following problem, which should be quite easy to solve with Skyfield, yet I am a bit overwhelmed by all the different coordinate systems and the translation between them.

I have a Topos location on Earth and a Topos location of a LEO satellite. I am considering the line-of-sight between them. I want to determine the latitude and longitude of the position along this path, where it intersects a specific layer of the atmosphere.

An example would be the mesosphere and an existing dataset on its properties at around 100km that is given based on latitude and longitude. The intersection would allow me to better understand the interaction these properties have on the communication with the satellite.

I tried doing it with Skyfield directly, but only end up with an Apparent object that I cannot convert back to latitude, longitude on Earth. First, I trigonometrically determined the distance from Earth to the point, where the height of 100km is reached.

Then, I took the position on Earth and used the unchanged elevation, azimuth to keep the direction of the path and finally added the calculated distance to arrive at this position. I think I need to get a Geocentric object to use subpoint() in order to get the desired latitude, longitude of this location.

This is what I have so far:

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))

I have gone through the documentation multiple times, and stumbled upon frame_latlon(frame) for Generic ICRF objects, but am not sure how to further proceed.

Trying it completely trigonomatrically with the latitudes and longitudes didn't yield the desired results either.

Unfortunately I do not really have any validated results that could be used to solve this problem more easily. Imagining it again trigonometrically, it is obvious that an increase in altitude of the satellite position would move the lat, lon of the intersection closer to the position on Earth. Decreasing the altitude would then move this intersection closer to the satellite.

Answer 1

That is an interesting problem, which Skyfield’s API provides no easy way to ask about; if you could outline the larger problem that will be solved by knowing the intersection of the line-of-sight with a particular altitude, then it is possible that a routine addressing that problem could be written for future users tackling the same question.

In the meantime:

1. To get your script running I had to import Distance from api.
2. The name dis was not recognized, so I replaced it with distance_to_height, hoping that it was the name intended.
3. Calling ts.now() is giving you a slightly different date and time on each call. While the script runs so fast that it probably does not matter, I have for clarity pivoted to calling now() only once at the beginning of the script, which is also slightly faster than calling it repeatedly. (Actually in this case it’s much faster, because the rotation matrices only get computed once rather than having to be computed over again for each separate time object, but that’s a hidden detail that’s not easy to see.)
4. I suspect a problem with your geometry: the 100 / sin() maneuver would only work if the Earth were flat, I think? But maybe you are always dealing with nearly-overhead satellites and so the error is manageable? (Or maybe I am mis-imagining the geometry; feel free to provide a diagram if the math is in fact correct.)
5. For readability I’ll give the components of altaz() names rather than numbers.

With those tweaks in place, I think the answer is that you need to manually construct a Geocentric position by adding together the position of the observer and the relative vector you have created between the observer and the kind-of-100km point along the line of sight. Having to take a manual step like this suggests a possible area where Skyfield can improve. Here is how it looks in code:

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')

The result I see is:

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

And for the future, I have added some thoughts to the Skyfield TODO.rst file that together might move towards unlocking a more idiomatic way to perform this kind of calculation in the future — though I suspect that a few more steps even beyond these will be necessary:

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