# DiscussionImproved (Tesseral) Gravity Model

#### n72.75

Tutorial Publisher
Donator
I updated my branch with the latest commits, and thought it would be a good opportunity to a record a (relatively) short video of this in action.

GLS

#### n72.75

Tutorial Publisher
Donator
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.

#### Ajaja

##### Active member
I believe that the differences result from a difference in the moon's orientation between GMAT and Orbiter.
I don't think so. The difference in the moon's orientation for MJD 52006 between GMAT's (MOON_PA) and Orbiter ~0.001 rad.
Have you compared results with GMAT\data\gravity\luna\grgm900c.cof model?

Last edited:

#### n72.75

Tutorial Publisher
Donator
I don't think so. The difference in the moon's orientation for MJD 52006 between GMAT's (MOON_PE) and Orbiter ~0.001 rad.
Have you compared results with GMAT\data\gravity\luna\grgm900c.cof model?
The position error on the first timestep was ~2.3m which agrees with that factor, roughly.

I haven't tried this test with grgm900c yet, running these tests takes an hour or so a piece so producing a lot of test data is slow going.

What would probably be good to know is the other sources of error, and the maximum position offset I could expect from sources other than the gravity models.

#### jarmonik

##### Well-known member
Orbiter Contributor
Beta Tester
I am kinda surprised by the huge impact of this gravity model. Orbit that's initially circular fluctuates over 10km in altitude over mean radius, that's a lot. It's going to make lunar orbit navigation a lot more complicated. I guess that Olympus on Mars is also going to mess up with gravity.

#### n72.75

Tutorial Publisher
Donator
The effect is quite dramatic at low altitudes, which you can easily reach due to the lack of atmosphere.

#### n72.75

Tutorial Publisher
Donator
Some slight improvement:
I changed the moon mass (M = GM/G) in the moon.cfg to be identical to that reported in LP165, also I improved the lunar equatorial radius with a bit more precision.

This now results in a graph that is so nice, that it's hard to see how nice it really is (compare to above):

Unfortunately these two are still up to their usual tricks.

More models need to be tried...

#### jarmonik

##### Well-known member
Orbiter Contributor
Beta Tester
Since the drift rate is very linear it's not a force we looking for more likely a velocity offset. You mentioned something about first time step but unfortunately I don't know much about it. However, I have noticed an erratic velocity error in the Moon that's pretty large. That can be reduced by increasing the precision from the Moon.cfg if you haven't tried that yet then it might be good idea to try. Also a reference frame conversion could be the source of the problem. I am not familiar with GMAT but I would assume it using ICRF.
I once (10 years ago) made the Moon ephemeris module using DE405 which should also use ICRF, I'll look if I can find it.

#### jarmonik

##### Well-known member
Orbiter Contributor
Beta Tester
Found it. It's 11.7MB in size. Can you take a file of that size by e-mail ? Also you may need to recompile it to switch off approximated ecliptic conversion. DE405 files did not come with a license so, for further distribution the license should be verified from JPL.

#### n72.75

Tutorial Publisher
Donator
I did bump up the precision limit in moon.cfg so it now uses all of the ELP82 terms that Orbiter has.

GMAT, at least the public 2020 version defaults to DE405. That would be useful to compare. I will PM you my Gmail. Thank you.

I should also verify that my parser is working correctly and that it's not stripping of the negative signs after row 10000 or anything weird like that. It's based on scanf and Ive tested it quite thoroughly, so I doubt it.

#### Ajaja

##### Active member
I would also recommend SPICE module:
I've uploaded the latest ready-to-use build with all configs and SPICE kernels here:
For x64 version of Orbiter you need to rename spice_x64.dll to spice.dll.
And you can use SPICE+orbiter.bsp kernel in GMAT's SolarSystem parameters to get the same position in GMAT and Orbiter.
orbiter.bsp is a DE440-based kernel (the spkmerge file - Kernels\orbiter.mrg).

And if you need the same orientation of the Moon there is the Orbiter 's parameters for MJD 52006, they are calculated from MOON_PA reference frame that GMAT uses:
Code:
SidRotPeriod =  2359031.191366999
SidRotOffset =  5.127420079135963
Obliquity =  0.026288976972756034
LAN =  1.7553689545171214
LAN_MJD =  51544.5
PrecessionPeriod =  0
PrecessionObliquity =  0
PrecessionLAN =  0

#### Ajaja

##### Active member
And one more very important thing, the MJD epoch in Orbiter is TDBModJulian+29999.5 epoch in GMAT. Not UTC/TT/A1/etc.

#### n72.75

Tutorial Publisher
Donator
And one more very important thing, the MJD epoch in Orbiter is TDBModJulian+29999.5 epoch in GMAT. Not UTC/TT/A1/etc.
This is a very good thing to know. Unfortunately, this does not magically fix everything like I was hoping when I first read it.

I've tried a number of tests with the SPICE modules. They drop in and run very nicely into my OpenOrbiter build.
Using SPICE has improved my initial position error so it is now [0.38335, -2.8685, -0.48492] mm

Unfortunately there does still seem to be a velocity offset somewhere, and it does seem to be somewhat worse, something like 8mm/sec vs about 5mm/sec. I'm sure that the SPICE modules are vastly more accurate than ELP82. The agreement I am trying to get between GMAT and Orbiter is more to demonstrate that the implementation of Pines algorithm, doesn't have any fundamental errors. I am very hesitant to pull request something that may have an underlying error in it, because the bugs it could cause down the road for people could be all but impossible to solve.

Interestingly enough, If you completely disable non-spherical gravity in both GMAT and Orbiter, the results are strikingly similar.

I am fairly confident that the drift in position has nothing to do with my gravity model. Things like sign errors on position and acceleration cause rapid accumulation of hundreds of kilometers of error, likewise with reference frame errors (mul/tmul confusion).

I don't have a good grasp of how error control works in either Orbiter or GMAT, but that's probably my next thing to look at.

#### n72.75

Tutorial Publisher
Donator
Subtracting the error in the spherical data from the non-spherical data looks like this:

#### n72.75

Tutorial Publisher
Donator
Vesta has a cute little 20x20 SHADR model. Very fun to fly in.

Replies
0
Views
122
Replies
0
Views
198
Replies
11
Views
714
Replies
11
Views
1K
Replies
180
Views
5K