Change the orbital analyser to use the Trajectory interface

[pleroy]

(There is still one TODO.)

#3493

🍐

Run on (4 X 3311 MHz CPU s)
CPU Caches:
  L1 Data 32 KiB (x4)
  L1 Instruction 32 KiB (x4)
  L2 Unified 256 KiB (x4)
  L3 Unified 6144 KiB (x1)
----------------------------------------------------------------------------------------------------
Benchmark                                                          Time             CPU   Iterations
----------------------------------------------------------------------------------------------------
OrbitalElementsBenchmark/ComputeOrbitalElementsEquatorial  579610000 ns    578125000 ns            1

[eggrobin]

Notes from various experiments on this test:

• increasing the number of steps per period has no effect;
• if downsampling is removed, the results are consistent with the old expectations, both with the old and the new OrbitalElements logic;
• the downsampling tolerance is 10 m.

[pleroy]

With a tolerance of 10 m, the trajectory has 149 776 points.

With a tolerance of 1 m, the trajectory has 259 023 points, but the semimajor axis range is only 14 909.5 km .. 14 910.7 km. So we should just lower the tolerance for the orbit analysers.

[eggrobin]

Further investigation: at 10 cm the range becomes 14 909.9 km .. 14 910.3 km, which seems like a relevant difference given the size of the range under consideration. Lowering the tolerance further does not change the range to a precision of 100 m; in particular that range is the same in the absence of downsampling.

There is a separate issue, which is that the estimated range is highly sensitive to the value of 24 points per period (with no downsampling, 96 points gives a range of 14 909.7 km .. 14 911.0 km, and 120 points a range of 14 910.5 km .. 14 911.0 km. This is because, contrary to the apsides, we make no attempt to interpolate when estimating the extrema.

This is a problem, but it is a problem we had already (we got our arbitrary samples from the integration of motion followed by downsampling instead of getting them from the averaging of elements, but there was no attempt at interpolation either). We should now have the tools to tackle it properly, since we are already doing most of the work needed to compute the derivatives of the mean elements; but this is a matter for another pull request.

[eggrobin]

it is a problem we had already

Hm, no, this does not make sense. We would not be getting maxima larger than the correct one nor minima smaller than the correct one just because of poor sampling. Is the integrator not converged enough?

[eggrobin]

Is the integrator not converged enough?

Indeed; the integrators converge only to order 1 or 2 (at least at the precisions that are relevant here), because the of the discontinuities in the higher derivatives of the downsampled trajectory.

Aside: it is clear that here we want high enough precision that a dozen points per period is not enough, and for an elliptic orbit that means we benefit from a variable stepsize.

Possibilities: Heun–Euler (but that is unsatisfactory, it is mostly cutting our losses), or lowering the downsampling tolerance (maybe to 1 mm?) until Dormand–Prince converges properly.

[pleroy]

Call this t_mid?