NaN derivatives in Keplerian propagation with Cartesian coordinates

In FieldCartesianOrbit, shiftPVElliptic causes NaN on the derivatives when used with Gradient (although this probably happens as well with DerivativeStructure, etc.) and a circular orbit, in particular of the type position (X,0,0) and velocity (0, sqrt(mu/X),0). This is not a singularity of the propagation itself (as there is no problem if converting to equinoctial, propagating and converting back), but rather of the particular method used. The NaNs appear at the call to arctan2