Based on @bgodard's comments on #45 it seems we still do not have an optimal design. The current implementation probably works for anybody not doing interplanetary flight dynamics (modulo bugs). Addressing this topic will probably entail another redesign and this time I want to get it right. Thus, I would like to collect some user stories and test cases here.

From #45

julia> t0=TAIEpoch(2000,1,1,12,00,32.0)
2000-01-01T12:00:32.000 TAI

julia> t1=UTCEpoch(2000,1,1,12,00,32.0)
2000-01-01T12:00:32.000 UTC

julia> t2=UTCEpoch(0.0,t0)
2000-01-01T12:00:32.000 UTC

julia> t2-t1
0.0 seconds

julia> t2==t1
false

julia> t2>t1
false

julia> t2<t1
false

julia> t2,t1
(2000-01-01T12:00:32.000 UTC, 2000-01-01T12:00:32.000 UTC)

I have also decided to restrict the input to TAIEpochs when overriding S-TAI.

This is not convenient at all because transformation between TDB and spacecraft proper time in a deep space trajectory has no meaningful relation to TAI. User would have to invent a fake intermediate TAI transformation.

julia> t2=TDBEpoch(2000,1,1)
2000-01-01T00:00:00.000 TDB

julia> t1=TDBEpoch(2000,1,1,12)
2000-01-01T12:00:00.000 TDB

julia> t2-t1
-43199.99998551589 seconds

But I would like to get the number of TDB seconds not TAI (TAI difference is not very useful for solar system dynamics though this is a good approximation for most applications).

Maybe we could have difference of Epoch in same timescale return seconds in that timescale.
For mixed timescales Epoch difference, raise an Error or convert the second epoch argument to the timescale of the first epoch argument (could be explicitly enabled by "using UnsafeEpochConversions")

At the moment I see it as a big problem that all timescales are directly related to TAI. I would like to have possibility for epochs which do not have any valid default conversion to TAI (timescale graph is disconnected), but possibly to other custom timescales.