title: Path between two Topos() locations: determine latitude and longitude where a given altitude is reached
body: <p>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.</p>
  <p>I have a <code>Topos</code> location on Earth and a <code>Topos</code> 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.</p>
  <p>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.</p>
  <p>I tried doing it with Skyfield directly, but only end up with an <code>Apparent</code> 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.</p>
  <p>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 <code>Geocentric</code> object to use <code>subpoint()</code> in order to get the desired latitude, longitude of this location.</p>
  <p>This is what I have so far:</p>
  <pre><code>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))
  </code></pre>
  <p>I have gone through the documentation multiple times, and stumbled upon <code>frame_latlon(frame)</code> for Generic ICRF objects, but am not sure how to further proceed.</p>
  <p>Trying it completely trigonomatrically with the latitudes and longitudes didn't yield the desired results either.</p>
  <p>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.</p>