Look at that! Back in 2015 someone contributed an example routine that
I just re-wrote without remembering it was there. I guess I should grep
for “new” features before sitting down to re-implement them.
The routine would die because `combine()` is a class method, not a
module method; and because it needed a timezone. A test is here added
to keep the feature fixed. Also, the timelib module used `dt` for the
`datetime` module but also sometimes as a convenient local name for a
`datetime`, so let’s rename the module to `dt_module` for readability
and consistency. Inspired by the attempted fix in #410.
This eliminates repeated deprecation warnings when rendering the docs:
```
MatplotlibDeprecationWarning: Adding an axes using the same arguments as
a previous axes currently reuses the earlier instance. In a future
version, a new instance will always be created and returned. Meanwhile,
this warning can be suppressed, and the future behavior ensured, by
passing a unique label to each axes instance.
```
Because two components of the position vector are used to finish
building the velocity vector, it was a mistake to rotate the position
vector using the M matrix until after the velocity vector has also been
built. This increases by almost two magnitudes the agreement between
Earth satellite reported velocities and the actual difference between
their successive positions.
This eliminates one source of noise as I dive into why Skyfield
satellite velocities do not hew as closely to actual satellite motion as
they do when the underlying sgp4 library is used alone.
It turns out I was somehow being hamfisted in how I built the Julian
date to pass to `TEME_to_ITRF()` in the test. By building the Julian
date manually, without going through a `Time` object, I can produce far
more exactly the quantities specified in AIAA-2006-6753 Appendix C.