Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change the orbital analyser to use the Trajectory interface #3519

Merged
merged 65 commits into from
Feb 12, 2023
Merged

Change the orbital analyser to use the Trajectory interface #3519

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
65 commits
Select commit Hold shift + click to select a range
5d2e8cd
Start to use Trajectory, but unwinding makes our life harder.
pleroy Jan 1, 2023
600bbf0
A promising attempt
eggrobin Jan 1, 2023
af60f79
the integrator is confused
eggrobin Jan 2, 2023
f0bfe5f
merge
eggrobin Jan 2, 2023
95a1d7b
unwinding
eggrobin Jan 2, 2023
c95b242
unwinding that works
eggrobin Jan 2, 2023
437f173
more traces
eggrobin Jan 2, 2023
788ca62
back to fixed step
eggrobin Jan 2, 2023
0e66b9a
bad unwinding
eggrobin Jan 2, 2023
8400e7e
An attempt at Clenshaw-Curtis that does not work too badly
eggrobin Jan 2, 2023
c3501d4
do we really need dense output?
eggrobin Jan 2, 2023
5414a4d
Merge remote-tracking branch 'la-vache/master' into better-integration
eggrobin Jan 12, 2023
ca81bbd
use RK4
eggrobin Jan 12, 2023
d19fbd3
a bug
eggrobin Jan 15, 2023
0d6dbd6
Trying to throw more automated Clenshaw Curtis at the problem
eggrobin Jan 16, 2023
b9ca286
do interpolation correctly
eggrobin Jan 16, 2023
8ac094d
compute osculating elements, and notice the noise
eggrobin Jan 19, 2023
50959b2
Merge branch 'better-integration' of https://github.com/eggrobin/Prin…
pleroy Jan 28, 2023
5e8ed38
ELMS.
pleroy Jan 28, 2023
e5d28e6
Merge.
pleroy Jan 28, 2023
2de7b6d
Merge branch 'Orbit' into Benchmarque
pleroy Jan 28, 2023
9bc51ed
AB order 6.
pleroy Jan 28, 2023
44bffc8
Merge branch 'Benchmark' into Benchmarque
pleroy Jan 28, 2023
afcab79
remove logger and don’t densely sample outside of tests
eggrobin Jan 28, 2023
024d755
Avoid allocations in the integrator.
pleroy Jan 29, 2023
ad1eebd
Avoid allocation in the ODE's equation.
pleroy Jan 29, 2023
b3ffd31
Unused code.
pleroy Jan 29, 2023
c9c0dcd
Move instead of copy.
pleroy Jan 29, 2023
4f40081
Merge.
pleroy Jan 29, 2023
5ca5aa3
Adjust for non-vector ODEs.
pleroy Jan 29, 2023
6286a1f
Je n'en sais rien.
pleroy Jan 30, 2023
aaad34a
Let's make a pause.
pleroy Jan 30, 2023
536f97b
A benchmark with a more inclined trajectory.
pleroy Jan 30, 2023
38970e6
Cleaned up the Augean stables.
pleroy Jan 30, 2023
c100c60
It kind-of works.
pleroy Jan 30, 2023
0a31db0
All tests passing.
pleroy Jan 30, 2023
58440f8
Cleanup.
pleroy Jan 30, 2023
d866566
Bring the projects back alive.
pleroy Jan 30, 2023
6fc31b8
Fix a compilation error.
pleroy Jan 30, 2023
2342674
Minor fixes.
pleroy Jan 30, 2023
04b2c47
Tolerances.
pleroy Jan 31, 2023
8f3720c
Lint.
pleroy Jan 31, 2023
1c290e3
Parameterization of Clenshaw-Curtis.
pleroy Jan 31, 2023
47a87f5
Specific downsampling tolerance for the orbit analyser.
pleroy Feb 7, 2023
bbe0517
Smaller tolerance.
pleroy Feb 7, 2023
413773c
Fix compilation error.
pleroy Feb 7, 2023
89129e7
Merge branch 'Benchmarque' of https://github.com/pleroy/Principia int…
pleroy Feb 8, 2023
35612b0
Heun
eggrobin Feb 8, 2023
cef4f0b
remove bad condition
eggrobin Feb 8, 2023
5629624
Compile out the logging.
pleroy Feb 8, 2023
86167ae
Mathematica logging.
pleroy Feb 10, 2023
ec5d1c5
Moal traces.
pleroy Feb 10, 2023
df94c90
Merge.
pleroy Feb 10, 2023
5127615
Code cleanup and better tolerances.
pleroy Feb 10, 2023
778f2eb
Bring back the tests.
pleroy Feb 10, 2023
a547281
Fix the tolerances.
pleroy Feb 10, 2023
9b6c6dd
Fixed-step integrator in test.
pleroy Feb 10, 2023
778d1a8
Minor cleanup.
pleroy Feb 11, 2023
c10ba43
More funky tolerance.
pleroy Feb 11, 2023
836ad9e
Logging the errors.
pleroy Feb 11, 2023
34d3ebc
Combined.
pleroy Feb 11, 2023
98cbeb8
Remove bad tests.
pleroy Feb 12, 2023
6302dc3
Benchmark downsampling tolerance.
pleroy Feb 12, 2023
b30a10b
After egg+phl's review.
pleroy Feb 12, 2023
3d4be26
Compilation error.
pleroy Feb 12, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 11 additions & 9 deletions astronomy/orbital_elements_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,7 @@ using testing_utilities::IsNear;
using testing_utilities::IsOk;
using testing_utilities::RelativeError;
using testing_utilities::operator""_;
using ::testing::AnyOf;
using ::testing::Lt;

class OrbitalElementsTest : public ::testing::Test {
Expand Down Expand Up @@ -187,16 +188,16 @@ TEST_F(OrbitalElementsTest, KeplerOrbit) {
// Mean element values.
EXPECT_THAT(elements.mean_semimajor_axis_interval().midpoint(),
AbsoluteErrorFrom(*initial_osculating.semimajor_axis,
Lt(410 * Micro(Metre))));
Lt(4.2 * Milli(Metre))));
EXPECT_THAT(elements.mean_eccentricity_interval().midpoint(),
AbsoluteErrorFrom(*initial_osculating.eccentricity,
Lt(1.9e-10)));
Lt(3.3e-10)));
EXPECT_THAT(elements.mean_inclination_interval().midpoint(),
AbsoluteErrorFrom(initial_osculating.inclination,
Lt(9.5 * Micro(ArcSecond))));
EXPECT_THAT(elements.mean_longitude_of_ascending_node_interval().midpoint(),
AbsoluteErrorFrom(initial_osculating.longitude_of_ascending_node,
Lt(106 * ArcSecond)));
Lt(116 * ArcSecond)));
EXPECT_THAT(elements.mean_argument_of_periapsis_interval().midpoint(),
AbsoluteErrorFrom(*initial_osculating.argument_of_periapsis,
Lt(202 * ArcSecond)));
Expand Down Expand Up @@ -299,7 +300,7 @@ TEST_F(OrbitalElementsTest, J2Perturbation) {
EXPECT_THAT(theoretical_ωʹ, IsNear(14_(1) * Degree / Day));

EXPECT_THAT(RelativeError(theoretical_Ωʹ, elements.nodal_precession()),
IsNear(0.0026_(1)));
AnyOf(IsNear(0.0026_(1)), IsNear(0.0029_(1))));

// Mean element values. Since Ω and ω precess rapidly, the midpoint of the
// range of values is of no interest.
Expand All @@ -314,21 +315,22 @@ TEST_F(OrbitalElementsTest, J2Perturbation) {

// Mean element stability: Ω and ω precess as expected, the other elements are
// stable.
EXPECT_THAT(elements.mean_semimajor_axis_interval().measure(),
IsNear(210_(1) * Milli(Metre)));
EXPECT_THAT(
elements.mean_semimajor_axis_interval().measure(),
AnyOf(IsNear(157_(1) * Milli(Metre)), IsNear(210_(1) * Milli(Metre))));
EXPECT_THAT(elements.mean_eccentricity_interval().measure(),
IsNear(3.7e-8_(1)));
AnyOf(IsNear(2.9e-8_(1)), IsNear(3.7e-8_(1))));
EXPECT_THAT(elements.mean_inclination_interval().measure(),
Lt(212 * Micro(ArcSecond)));
EXPECT_THAT(
RelativeError(
-theoretical_Ωʹ * mission_duration,
elements.mean_longitude_of_ascending_node_interval().measure()),
IsNear(0.0009_(1)));
AnyOf(IsNear(0.0009_(1)), IsNear(0.0023_(1))));
EXPECT_THAT(
RelativeError(theoretical_ωʹ * mission_duration,
elements.mean_argument_of_periapsis_interval().measure()),
IsNear(0.0018_(1)));
AnyOf(IsNear(0.0018_(1)), IsNear(0.0025_(1))));

mathematica::Logger logger(
SOLUTION_DIR / "mathematica" / "j2_perturbed_elements.generated.wl",
Expand Down
4 changes: 2 additions & 2 deletions astronomy/лидов_古在_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,9 +159,9 @@ TEST_F(Лидов古在Test, MercuryOrbiter) {
// is excellent: while the sun is nudging and deforming the orbit, it is not
// pumping energy into nor out of it.
EXPECT_THAT(elements.mean_semimajor_axis_interval().min,
IsNear(14'910.0_(1) * Kilo(Metre)));
IsNear(14'906.7_(1) * Kilo(Metre)));
EXPECT_THAT(elements.mean_semimajor_axis_interval().max,
IsNear(14'910.3_(1) * Kilo(Metre)));
IsNear(14'913.5_(1) * Kilo(Metre)));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member Author

@pleroy pleroy Feb 1, 2023 •
edited by eggrobin
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


// The integral c₁ is preserved quite well: we have an exchange between
// inclination and eccentricity.
Expand Down