As we near the one year anniversary of my original post in this thread, I figure it's time for some kind of update.
In the interest of making this feature worthwhile and something that improves the overall accuracy of the simulation, and not just something that messes up your orbit just to make things more difficult, I've been working on quantifying errors, and I know how much everyone likes my graphs.
My tool of choice, mostly because of ease and familiarity is GMAT and Octave(MATLAB).
The test goes like this:
- Set up an orbit in Orbiter and GMAT.
- Propagate for a while, collecting position data.
- Compare results
- Decide what an acceptable amount of error is.
The orbit chosen for this test:
- MJD 52006.76972
- 1800 km semimajor axis
- 1E-7 eccentricity
- 45° inclination
- 0° longitude of ascending node
- 0° true anomaly at epoch
This orbit is defined by this state vector in moon rotating equatorial;
MJD 52006.76972 (22007.26558510206 GMAT MJD)
X 1800000.018000013
Y -2.359001882723533e-9
X 1.77635683940025e-9
Vx 3.654845523448635e-13
Vy 1167.000010008068
Vz 1167.000010008068
Both test cases use the LP165 spherical harmonic model with all 165x165 terms, a 5 second timestep and 8th order Runge-Kutta propagator.
The results I'm getting are certainly not in 100% agreement with GMAT, but they are very close. To the best of my ability, I believe that the differences result from a difference in the moon's orientation between GMAT and Orbiter.
Shown below is a plot of altitude above the mean lunar radius, comparing my Pines implementation to the implementation found in GMAT.
The general trend of changing eccentricity, periapsis and apoapsis height over time are in agreement between the two software packages.
To quantify the error in position we can compare absolute positional error over time between the two packages. This is about two orders of magnitude larger than I would like, but it's likely not the fault of the gravity model itself.
Looking at the radial component of position error only, we see it is much more mild and approximately an order of magnitude less than the total error, indicating that the majority of the errors lie along the longitude-latitude directions.
~~~~
It is easy to be a bit dis-heartened by this error, but when we compare GMAT-LP165 against Orbiter's moon, with spherical gravity (as is the case in Orbiter 2016 and prior), we see a substantial difference in the error evolution over time.
With a spherical lunar gravitational potential field this state vector results results in a nearly circular orbit, the only perturbations upon which are the Earth and the Sun.
The way orbiter calculates its first timestep actually results in a rather substantial initial offset, which initially settles down, but then climbs again past 640km difference by the conclusion of the test, with no signs of slowing.