We should really compute the collisions in the C++. How hard can it be to intersect a sphere and a segment?

the easiest way to find collision point - use bisect.

Also, planets are not exact spheres. They have mountains and valleys. So using bisect looks like the only solution

Planets are spheres in Principia, just say no to ridiculously expensive height maps. If your flight plan entails flying through the Grand Canyon, think again because you'll be sorry as soon as you enter the atmosphere.

Notes to self:

1. Trajectories are cubic, so the intersection involves solving a cubic equation. That's hard to do in a way that's properly conditioned. There are papers by Kahan and Blinn on that topic.
2. We do use Bisect to find ascending/descending nodes. That's stupid. We should either use Brent or solve the cubic.

In the end, we decided to compute the impact properly, taking the terrain into account (and not assuming that planets are spheres).

For the record, it's not practical to find the the first collision using bisection. Bisection will easily find one collision, but not necessarily the first. In addition, in order to use bisection in the first place you need a point underground, and there is no guarantee that the periapsis is underground even if there is a collision: the periapsis could be in a crater and the collision on the wall of that crater. We used such an approach in Jordan and Julia, but it resulted in collisions being drawn at semi-random places or missed entirely.

In the end we found that building interpolants to the height above the terrain using Чебышёв polynomials, with suitable subdivisions, is the best approach. We then compute all zeros of these polynomials as eigenvalues of the companion matrix using the Schur decomposition, thereby ensuring that we find all the zeros regardless of the shape of the terrain. See the above pull requests and [Boy06] and [Boy13] in the bibliography.