# DiscussionImproved (Tesseral) Gravity Model

#### n72.75

##### Move slow and break things
Orbiter Contributor
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.

#### n72.75

##### Move slow and break things
Orbiter Contributor
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

##### Move slow and break things
Orbiter Contributor
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

##### Move slow and break things
Orbiter Contributor
Tutorial Publisher
Donator
The effect is quite dramatic at low altitudes, which you can easily reach due to the lack of atmosphere.

#### n72.75

##### Move slow and break things
Orbiter Contributor
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

##### Move slow and break things
Orbiter Contributor
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

##### Move slow and break things
Orbiter Contributor
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

##### Move slow and break things
Orbiter Contributor
Tutorial Publisher
Donator
Subtracting the error in the spherical data from the non-spherical data looks like this:

#### n72.75

##### Move slow and break things
Orbiter Contributor
Tutorial Publisher
Donator
Vesta has a cute little 20x20 SHADR model. Very fun to fly in.

#### n72.75

##### Move slow and break things
Orbiter Contributor
Tutorial Publisher
Donator
Good news everyone!

This has been merged into the main branch.

Go have fun with the latest build, don't forget to enable non-spherical gravity.

Because the perturbation force calculated by this model is done in a very similar way (same spot in the code), as the old zonal-only method, Orbiter's substeps and propagators work with it exactly the same as thy did before.

By default, all models (for planets/moons that have them) are set to 10x10 to minimize the impact on performance. This can be changed in the celestial body config file, along with the gravity model itself (as long as it's SHARD-ascii format).

Note:
• Increasing the GravCoeffCutoff value will impact performance with a time complexity O(n²)
• Increasing the number of vessels in a scenario is still O( n )...buuuuut if you want a huge gravity model with hundreds of terms, you may need to have less vessels, or vice versa. I tested a single Deltaglider in orbit around the moon with a 1200x1200 model, and it pulled my framerate down to about 1-2fps; its theoretically possible but there is really no need to ever go that high--propagator error will be higher than the millimeter precision you'd be after with that many terms.
• You don't need to download any gravity models to make this work; they are included with Orbiter. A table of bodies that include tesseral gravity models is show below
• The old non-spherical gravity is still there as a fallback for planets that only have (and really only need) the "Jcoeff", zonal terms.

 Body Gravity Model File Maximum Cutoff Mercury jgmess_160a_sha.tab 160 Venus mod_shgj120p.a01 120 Moon jgl165p1.sha 165 Mars jgmro_120f_sha.tab 120 Vesta JGDWN_VES20H_SHA.TAB 20

The source for all of these models is https://pds-geosciences.wustl.edu/ The Magellan Gravity Model was modified slightly to work with the header parser for the other models. Earth isn't on the list yet; the others all have far bigger perturbations, and I need to convert one of the available Earth models into SHADR-ASCII format.

There is still an underlying difference between Orbiter and reality, but I think this is more related to ephemerides and rotational models (more on this in a future episode), and its on the order of meters, not hundreds of kilometers.

This has been an enormous undertaking, and of the work required to make this happen, I have only done a very small part of the heavy lifting compared with those whom I have cited in documentation.

Have fun!

#### Attachments

• 1684455307737.png
33.7 KB · Views: 5

#### misha.physics

##### Well-known member
Hi,
I just wanted to clarify, if this non-spherical gravity model is included "as default" in the latest Orbiter build on github both for x86 and x64, and replaces original gravity model? And this model won't work when "non-spherical gravity sources" is disabled in "physics settings" tab, isn't it?

#### n72.75

##### Move slow and break things
Orbiter Contributor
Tutorial Publisher
Donator
Hi,
I just wanted to clarify, if this non-spherical gravity model is included "as default" in the latest Orbiter build on github both for x86 and x64, and replaces original gravity model? And this model won't work when "non-spherical gravity sources" is disabled in "physics settings" tab, isn't it?
You are somewhat correct.

The new model is included by default. It's not an addon, and it's tightly integrated into the Orbiter code the same way the old one is/was.

If (for instance for Jupiter) you only specify the jcoeff values, then Orbiter will use the old model, so it's fully backward compatible. But to be clear. The two don't work together: Orbiter either uses the old Zonal model, if jcoeff values are specified, and no gravity model file is loaded, or it loaded the new Tesseral model. This set on a body by body basis in the planet/moon config file, e.g. moon.cfg

Disabling non-spherical gravity disables all non-spherical gravity just like it used to.

#### Ajaja

##### Active member
Excellent job!
I did some test and compared the new improved gravity model in OpenOrbiter with GMAT's models (even made a script for converting *.tab files to GMAT's *.cof files).
Orbiter and GMAT are very close now.
An orbit around Mars for example:

But once, only once, I got a strange issue. And wasn't able to reproduce this:

I don't know what the hell it was. I made tests in the console mode, no fuel in the ship (the script and FlightData.log in attachments), so it wasn't some misclick. I used jgmro_120f_sha.tab model.
Next test - all ok.
Any idea what it might be?

#### Attachments

• m1.scn
361 bytes · Views: 0
• FlightData.log
72.4 KB · Views: 2
GLS

#### kuddel

##### Donator
Donator
Wild guess: Maybe a clock/timing issue? Somthing like "System Clock got adjusted by NTP", so the Operating-System skipped (or rolled back) some time.
That might have caused calculations in Orbiter to have the same effect as acceleration.
...wild guess however.

Replies
25
Views
1K
Replies
7
Views
396
Replies
0
Views
335
Replies
0
Views
428
Replies
11
Views
1K