Fast error-controlling MOID computation for confocal elliptic orbits

We present an algorithm to compute the minimum orbital intersection distance (MOID), or global minimum of the distance between the points lying on two Keplerian ellipses. This is achieved by finding all stationary points of the distance function, based on solving an algebraic polynomial equation of $16$th degree. The algorithm tracks numerical errors appearing on the way, and treats carefully nearly degenerate cases, including practical cases with almost circular and almost coplanar orbits. Benchmarks confirm its high numeric reliability and accuracy, and that regardless of its error--controlling overheads, this algorithm pretends to be one of the fastest MOID computation methods available to date, so it may be useful in processing large catalogs.