As the time between each line is really small (milliseconds), it is
better to propagate first and use an ephemeris later.

The low point at tile minimum elevation was computed with a tiny
positive error, meaning the point really was above min elevation (at
micrometer level). The line segment between entry and exit stopped just
before traversing the Digital Elevation Model, and an error was
triggered.

The solution was to take some margin when computing the segment
endpoints at entry and exit: we now start above max elevation and end
below min elevation so the intersection should be really on the line. 

There were some small errors (at centimeter level) as result points may
be slightly out of line of sight. These errors seemed to be due to the
final linear line-of-sight model that is used at pixel level.

An iterative correction step is performed to ensure the point is really
on the line-of-sight.

The residual error is now at nanometer level.

The conversion is of course accurate only in the neighborhood of the
reference point, as a straight line in Cartesian is not a straight line
in geodetic coordinates. What is converted is the initial direction so
the two curves are tangent at the reference point.

Near Earth surface, the two curves remain within one millimeter of each
other after about 100 m travel. 

The cases handled include linear cases (when Digital Elevation Model has
zero curvature) and constant cases (no curvature and line-of-sight
parallel to tile, both in or out-of-tile).

The case where the line-of-sight enters the Digital Elevation Model in
one tile on top and exit it in another tile on bottom is not tested yet.

The test does not succeed yet, an infinite loop is triggered.

These methods are essential to identify where to split the line-of-sight
in the Duvenhage algorithm, when going from one level in the min/max
kd-tree to the next level.