Max Degree / Order in GravityModel

Help
Anonymous
2013-03-14
2013-03-24
• Anonymous - 2013-03-14

Hi,

I want to use the GravityModel class to compute accelerations affecting satellites due to the gravitational potential. I wonder if it's possible to define the degree and order up to which I want the field to be calculated during the runtime. As far as I understand the code, it is always computing the gravitational field to the highest possible level?

• Charles Karney - 2013-03-14

At present a full sum is always done. You could construct model
files with a subset of the coefficients. Alternatively modifying
the code to allow the max degree and order to be specified in the
constructor would be possible. Actually, it looks really easy, just
modify the two calls to the SphericalHarmonic constructor in the
constructor for GravityModel and substitute smaller numbers for the
4th and 5th arguments. This isn't perhaps exactly what you want because
you would need to instantiate a new GravityModel for each set of
[max-degree, max-order] you want. But it might be enough to get you
started.

• Freddy Kay - 2013-03-15

Hi Charles,

Thank you for your help! Though I actually wanted to avoid modifiyng the code, I tried you suggestion anyway and it works just fine.
As far as I can tell there are only two calls to SphericalHarmonic, just as you said, or did I miss one?
And another question: have you ever crosschecked the results of your Code? Before I tried your library, I already used another implementation which can be found here (see GEOGRAV). The results I get differ by aproximately 0.1%. Do you have any suggestions as to why this is the case?

FYI:
Method: GravityModel::V();
ECI-Coords used [m]:
X = -5522333.38232309840
Y = 4238164.28094564470
Z = 1208232.17500043990

Results:
gx = 6.2485000069713497
gy = -4.7954346178341929
gz = -1.3708162469637046

Results by GEOGRAV (though here they use ECEF Coordinates):
ax = 6.2481691347188901
ay = -4.7954745597432824
az = -1.3707634134685798

• Charles Karney - 2013-03-15

I guess that you are using the egm96 gravity model??

I carried out various checks on the Gravity class: specifically
verifying that the spherical harmonic sums are carried out correctly
(without underflow or overflow) and cross-checking the results it gives
for the geoid height against the tabulated geoid heights given by NGA.

As to what can cause the discrepancy of 0.1% between GeographicLib and
another software package... Possibilities are:

• accumulated round-off errors (but this is unlikely with egm96). I ran
the test point you gave comparing a calculation with double precision
(53 bit fractions) and long doubles (64 bit fractions). Calling
V(-5522333.3823230984L, 4238164.2809456447L, 1208232.1750004399L,
GX, GY, GZ);
gives for GX, GY, GZ
53-bit 6.2485000069713497 -4.7954346178341929 -1.3708162469637046
64-bit 6.2485000069713498 -4.7954346178341933 -1.3708162469637048
So GeographicLib is calculating the sum accurately. You can try
reducing the degree and order for both packages to check whether
round-off is a problem with GEOGRAV.

• your mysterious distinction between ECI coordinates (which I've never
heard of) and ECEF coordinates (which I would call geocentric).
You're on your own on this one.

• different gravity models. I assert that GeographicLib does things the
NGA way; but let me know if this is not so. GEOGRAV uses a model that
is subtly different (different values for a, GM, C21, S21. The
differences don't quite look big enough to explain the 0.1%
discrepancy. But you can try fiddling with these constants in one or
other of the two packages to check.

You're quite right to worry about a discrepancy of this magnitude. In
your shoes I would want to isolate the reason for the discrepancy.

• Freddy Kay - 2013-03-20

Hey Charles,

I tried to figure out where the discrepency is coming from. I used STK to calculate the accelerations and fact is I get results which are in compliance with the results in GEOGRAV and show more or less the same discrepency which I mentioned above. As you already checked the sum calculation, I suspect there is a problem with the reference frames.

ECI (Earth centered inertial) and ECEF (Earth centered Earth Fixed) are both geocentric coordinates. THe ECI frame refers to a system where the z-Axis is pointing towards the northpole/earth rotation axis and the x Axis is pointing towards the vernal equinox, most commonly it is ecpressed in terms of cartesian coordinates but you can also use spherical coordinates (Right ascension, declination and distance). The ECEF frame on the other hand is a system where the z-Axis is again directed towards the northpole but the x Axis is pointing towards the Greenwich meridian, it is therefore rotating in coincidence with the Earth's rotation It is most commonly expressed in terms of longitude, latitude and height above ground/distance from origin.

Whats weird when I use your library is, that when I use coordinates in the ECI frame, the results I get are of the same magnitude. If I use an ECEF frame I cannot reproduce the results I get with your library. As an ECI frame would need some kind of information on the timeframe in order to evaluate the geopotential of a rotating body, it seems weird that the results of GeographicLib are in the same region as the results I get with Geograv and STK (which include a time dependent conversion).

I already started to go through your code to find where and how you are converting coordinates in GravityModel but figured it would be easier to ask right away. So, where are you converting coordinates in GravityModel and are you converting them at all? And more important which Reference System are you using? Does it account for the true centre of gravity or is it using the geometric centre of earth? Are you maybe using a Coordinate System that is refering to an Epoch?

I hope you can tell me where I might have misunderstood your library or confused the coordinates or tell me if I am completely off track.

Thanks a lot in advance!

• Charles Karney - 2013-03-24

The coordinates systems that GeographicLib uses are all tied to the
earth. The only conversions done are simple ones, geocentric to
geographic and the like. I believe I am faithfully implementing the
EGM96 and EGM2008 gravity models published by NGA. So you should
consult the NGA documentation for an explanation of the model:

I think you'll find that the results GeographicLib gives accurately
matches the results from the NGA programs, e.g., f477.f or
hsynth_WGS84.f.

Incidentally there's no time dependence incorporated into these models.

Anonymous