Numerical stability and angles
The codebase has several angle computations whose numerical stability has not been investigated.
In particular, https://github.com/mockingbirdnest/Principia/blob/69cf1a4/geometry/r3_element_body.hpp#L108
computes the latitude as arc sin(z / sqrt(x^2 + y^2 + z^2)), which is likely to cause trouble because of arc sin's singularity. arc tg(z / sqrt(x^2 + y^2)), computed with the 2-parameter arctangent, should be stabler, although a detailed analysis would be desirable.
In addition, https://github.com/mockingbirdnest/Principia/blob/69cf1a4/physics/ephemeris_test.cpp#L710 computes an angle between two vector as the 2-parameter arctangent of the norm of the cross product.
While this does not present the singularity-induced instability of the arccosine method, this is not the method recommended by Kahan in https://www.cs.berkeley.edu/~wkahan/Mindless.pdf (p. 47).
In any case, the angle of two vectors is a useful enough operation that it deserves a dedicated function in geometry: astronomical tests, e.g., @pdn4kd's eclipses, will need this.