Approximating mean sea level / geoid radius with WGS-84 spheroid/ellipsoid

[simfinite]

Hey all,
I am using JSBSim together with the ArduPlane framework to simulate fixed-wing drones. Looking at JSBSim CSV logs, I noticed that the ECEF position coordinates were way off from what I expected. The aircraft was initialized at ~9000m height above the WGS-84 ellipsoid (GPS altitude) when initialized with altitude AGL set to '0'.

Looking around the code, I found out that this is because the default mean sea level radius (in FGDefaultGroundCallback) is initialized to the major semi-axis (equatorial earth radius). This constant value is used as the default reference for altitude AGL and terrain elevation everywhere.. the earth's a sphere again.

Now, I realize that there may be other implementations of FGGroundCallback, possibly with complex geoid and terrain models. However, for my use-case, all I really want is to be able to place the aircraft on an airfield, take-off and produce some reasonable log files. The WGS-84 spheroid is a good reference for that, especially since GPS altitude values can be found everywhere.

Does anyone agree that the WGS-84 spheroid/ellipsoid would be a much better approximation of the geoid than a simple sphere? Any reasons against it? Anything I may have missed/misunderstood?

I would be happy to implement this, if we can agree that it would make sense and wouldn't break anything...

[simfinite]

After some further digging into the code, I realized that the spherical earth (or spherical mean sea level) assumption is quite deeply baked into the current implementation (e.g. with concepts such as a "geocentric terrain level radius" which is useful only under this assumption). I played around a bit to get a deeper understanding of what would have to be done to have both altitude and terrain level defined as a height above an ellipsoid instead of a radius of a sphere. I am currently stuck on getting the trim routine to work with my changes.
Anyway, before going further, I would appreciate any feedback on whether this has been adressed before and if any other parts of JSBSim / FlightGear rely on the spherical mean sea level assumptions.. Also, I would like to know how this is dealt with from the FlightGear side, where one can travel around a non-spherical earth right?

[bcoconni]

Hi @simfinite
The Earth mean radius is indeed used in several places in JSBSim for different purposes (off the top of my head: gravity computations, standard atmosphere model and altitude/location computations). This certainly needs some love to make sure that:

1. The mean radius value is the same everywhere which I am pretty sure is not the case at the moment.
2. The mean radius is only used where it makes sense which is not guaranteed at the moment.

Contributions are welcome and if you want to address this topic then be prepared to fight a beast 😉 because there is interdependence between these computations and the rest of the code. But it can be a very interesting subject nonetheless.

In the short term, the default version FGDefaultGroundCallback can be customized with the property position/terrain-elevation-asl-ft which defines the terrain elevation above sea level in feet. That should allow you to place the ground at the expected distance from the Earth center. Altitude above Ground Level (AGL) should be used in that case.

[bcoconni]

Does anyone agree that the WGS-84 spheroid/ellipsoid would be a much better approximation of the geoid than a simple sphere? Any reasons against it? Anything I may have missed/misunderstood?

Modeling the Earth as a sphere is certainly something that should be kept in JSBSim (as least as an option) since most flight dynamics textbooks are either assuming that the Earth is flat (which we are not modeling) or that the Earth is a sphere because the maths are still manageable which is not the case if you assume that the Earth is an ellipsoid or even worse a general geoid.

Having said that, nothing prevents us from adding an option where the Earth ground would be modeled as the WGS84 ellipsoid.

Also, I would like to know how this is dealt with from the FlightGear side, where one can travel around a non-spherical earth right?

FlightGear is using a specially designed derived class of FGGroundCallback which manages the very complex ground of the Earth. This class is used by JSBSim to interrogate FlightGear for the position of the ground without needing to know the details of the Earth shape (that's what OOP is all about after all).

However even when JSBSim is used within FlightGear, the sea level altitude is computed with the mean sea level radius that you mentioned above. I am not aware of anyone complaining that the location is inconsistent with GPS coordinates however (but this does not mean that the current code is correct, just that nobody noticed it before you did).

[bcoconni]

I dug into my personal archive and found much food for thoughts:

• This topic has already been raised by myself in 2013 😲 (see this mailing list thread - I completely forgot about it). The discussion is bringing many interesting points including some of those you raised. For some reasons, the issue made no further progress beyond this discussion.
• The bug report SF#80 is quite interesting in that it is raising a good point: the local frame might need some adaptations for the WGS84 ellipsoid. What this bug report had wrong is that the local frame should not be made tangent to the ellipsoid; instead the local Z frame should always be parallel to the gravity acceleration which it is currently not when using the J2 gravity.
• Yet another topic is: should the AGL (above ground level) altitude use geodetic coordinates ?
• Also the FGWaypoint class is using the Haversine formula to compute the heading and distance to another location on Earth. This would need to use another algorithm for the WGS84 ellipsoid.

I told you this is a beast to fight but this does not mean we cannot bring it down 😃

[simfinite]

Thank's for the feedback and digging up all these links. It's good to know that I seem to have hit a spot and not just failed to setup my simulation correctly. I agree it looks like a beast, but why not poke it a little bit to see what happens ;)

Before getting any deeper into the code and details of the issue, I would like to provide a user-side perspective of what I personally would like to have in JSBSim. As mentioned above, I use JSBSim in combination with the ArduPlane autopilot for simulating fixed-wing drones in the context of research projects. Sometimes, I need to export and work with the simulated trajectories in other tools/environments which is when the coordinate reference systems, earth models and altitude references really matter.

A simple example would be to export the trajectory to KML in order to visualize it in Google Earth. As explained above, if I would use current JSBSim's ECEF cartesian coordinates for this, the altitude would seemingly change with latitude due to the spherical earth in JSBSim and ellipsoidal earth in Google Earth. If I use lat/lon/alt coordinates from the spherical earth and project them on an ellipsoid it probably looks ok at first sight, but without figuring out the math of that projection and its implications, I don't really know what I'm doing anymore.

Here are the types of earth models I would like JSBSim to support along with the corresponding use cases:

• spherical earth: For simulations without requirements on accuracy and consistency of geodetic coordinates, e.g. simulating some dynamic effects, tuning flight controllers, simulating a final approach etc.. As I understand, this is implemented. The local terrain can be set with position/terrain-elevation-asl-ft and one can work with lat/lon and altitude AGL.

• ellipsoidal earth: For simulations which require accuracy and consistency of geodetic coordinates, e.g. flying over greater distances from one airport to another. Trying to do this with the currently implemented spherical earth results in the above mentioned problems.

• complex/geoid earth: For simulations which require accuracy and consistency of geodetic coordinates and complex terrain models. I understand that the complex terrain part is implemented in FlightGear's custom FGGroundCallback class. However, I am not yet convinced that the coordinate problems of the spherical earth are solved here, since you mention that "even when JSBSim is used within FlightGear, the sea level altitude is computed with the mean sea level radius".

From my current understanding we should distinguish between the two problems:

1. Providing ground level feedback (solved even for complex/generic case in FlightGear)
2. Being able to choose between a spherical and ellipsoidal mean sea level (with impact on atmosphere, gravity, altitude and distance calculations)

I would like to focus on the 2nd problem. Do you agree that this has not been solved for FlightGear yet? My current understanding is that FlightGear provides ground level ASL/MSL to JSBSim, which then maps this ground level onto its spherical mean sea level. If this is true, then the distance travelled in flights along a north/south axis on the spherical JSBSim earth should be slightly greater than expected for an ellipsoidal earth. Also, as mentioned above, the ECEF coordinates along this trajectory are not useful when applied in the context of an ellipsoidal earth model.

[bcoconni]

A simple example would be to export the trajectory to KML in order to visualize it in Google Earth.

There is already a similar feature in JSBSim. At the moment, it can only be generated via a script but it can certainly be moved to the input/output features of JSBSim.

jsbsim/src/input_output/FGScript.cpp

Lines 482 to 517 in 511869e

	// Print notification values after setting them
	if (thisEvent.Notify && !thisEvent.Notified) {
	if (thisEvent.NotifyKML) {
	cout << endl << "<Placemark>" << endl;
	cout << " <name> " << currentTime << " seconds" << " </name>" << endl;
	cout << " <description>" << endl;
	cout << " <![CDATA[" << endl;
	cout << " <b>" << thisEvent.Name << " (Event " << event_ctr << ")" << " executed at time: " << currentTime << "</b><br/>" << endl;
	} else {
	cout << endl << underon
	<< highint << thisEvent.Name << normint << underoff
	<< " (Event " << event_ctr << ")"
	<< " executed at time: " << highint << currentTime << normint << endl;
	}
	if (!thisEvent.Description.empty()) {
	cout << " " << thisEvent.Description << endl;
	}
	for (j=0; j<thisEvent.NotifyProperties.size();j++) {
	cout << " " << thisEvent.DisplayString[j] << " = " << thisEvent.NotifyProperties[j]->getDoubleValue();
	if (thisEvent.NotifyKML) cout << " <br/>";
	cout << endl;
	}
	if (thisEvent.NotifyKML) {
	cout << " ]]>" << endl;
	cout << " </description>" << endl;
	cout << " <Point>" << endl;
	cout << " <altitudeMode> absolute </altitudeMode>" << endl;
	cout << " <extrude> 1 </extrude>" << endl;
	cout << " <coordinates>" << FDMExec->GetPropagate()->GetLongitudeDeg()
	<< "," << FDMExec->GetPropagate()->GetGeodLatitudeDeg()
	<< "," << FDMExec->GetPropagate()->GetAltitudeASLmeters() << "</coordinates>" << endl;
	cout << " </Point>" << endl;
	cout << "</Placemark>" << endl;
	}
	cout << endl;
	thisEvent.Notified = true;

However, I am not yet convinced that the coordinate problems of the spherical earth are solved here, since you mention that "even when JSBSim is used within FlightGear, the sea level altitude is computed with the mean sea level radius".

I mentioned that off the top of my head but what I am sure of is that JSBSim does not depend by any means on FlightGear to perform its lat/lon/altitude computations. The FlightGear class FGGroundCallback is basically used to determine how far the ground is below the aircraft. For that purpose JSBSim supplies the aircraft geodetic coordinates and FlightGear sends back the height (which can be negative 😉) with respect to the ground.

I would like to focus on the 2nd problem. Do you agree that this has not been solved for FlightGear yet?

Well, somehow the FlightGear project has solved this problem and much beyond since it uses a complete and accurate model of the Earth geography based on "open source" cartography. However, even if JSBSim has tight links with FlightGear, it is an independent project meaning that we can have our own model of the Earth shape.
Of course, I would be reluctant to replicate a detailed Earth model such as FlightGear currently does but adding the WGS84 ellipsoid is OK to me (How could I pretend otherwise since I requested that myself a few years ago ? 😄 )

My current understanding is that FlightGear provides ground level ASL/MSL to JSBSim, which then maps this ground level onto its spherical mean sea level.

Honestly I do not know the details of the scenery implementation in FlightGear. What I know is that they use a completely separate project named TerraGear to produce and manage the data. I also know that TerraGear needs an insanely high amount of CPU/memory/disk space to run, up to the point that it has its own team devoted to the maintenance/development of the project. And since this is such a highly specialized and complex topic, I would certainly not want to duplicate their work.

If this is true, then the distance traveled in flights along a north/south axis on the spherical JSBSim earth should be slightly greater than expected for an ellipsoidal earth.

Well, the distance traveled only depends on the trajectory of the aircraft and not on the Earth shape. Having said that, the GNC (Guidance, Navigation and Control) that is generally implemented in JSBSim is indeed assuming that the Earth is a sphere meaning that the navigation results in a trajectory that might be slightly wrong as you mentioned. However JSBSim is so versatile that some people might have developed their own GNC that would account for the WGS84 ellipsoid (the Space Shuttle model comes to my mind as I know that its authors have come a long way towards fidelity).

Also, as mentioned above, the ECEF coordinates along this trajectory are not useful when applied in the context of an ellipsoidal earth model.

Indeed, the Space Shuttle model is the exception rather than the rule. Most users are using the standard features and I am convinced they would appreciate if JSBSim could save them the coding of the WGS84 computation in their model.

[simfinite]

Ok, thanks. I just wanted to make sure that I understand what has and what hasn't been solved in JSBSim and FlightGear.

the distance traveled only depends on the trajectory of the aircraft and not on the Earth shape.

Of course. But the trajectory of the aircraft is restricted to the atmosphere which again is bound to the surface (or actually the mean sea level) of the planet - spherical or ellipsoidal.

You mention implementation of GNC functions that rely on the spherical earth assumption. That probably includes the Haversine distance formula in FGWaypoint right? Any more examples?

[bcoconni]

@simfinite Sorry I lost track of this discussion but I have not forgotten that there is work to be done on this topic 😄

You mention implementation of GNC functions that rely on the spherical earth assumption. That probably includes the Haversine distance formula in FGWaypoint right?

You are correct this is what I was thinking about.

Any more examples?

Not off the top of my head but there must certainly exist some other examples.

BTW the work on this topic is also related to issue #201 since the crash described there is related to how JSBSim interacts with FGGroundCallback which is itself linked to the current discussion.

[bcoconni]

Resuming my work on this topic: it seems that we made a bad design decision years ago when making the distinction between geodetic and geocentric coordinates. Indeed there are 2 scenarios:

1. The Earth is assumed spherical in which case geocentric and geodetic coordinates are mathematically identical.
2. The Earth is assumed to be an ellipsoid in which case geocentric altitude is meaningless since the sea level radius of an ellipsoid is undefined.

In both cases the distinction is useless: a better decision would have been to use geodetic coordinates only and drop geocentric. Now we have to manage the situation where a user specifies a mix of geocentric latitude with geodetic altitude. The maths behind that scenario are not trivial and I am struggling to see the point of spending time implementing that (except for backward compatibility ?).

[seanmcleod]

Refreshing my knowledge on the difference between geodetic altitude and geocentric altitude, from - https://en.wikipedia.org/wiki/Geodetic_datum

geodetic altitude is defined as the height above the ellipsoid surface, normal to the ellipsoid; whereas geocentric altitude is defined as the height above the ellipsoid surface along a line to the center of the ellipsoid (the radius). When used without qualification, the term altitude refers to geodetic altitude; as is used in aviation. Geocentric altitude is typically used in orbital mechanics.

So something like this? GD-Alt (geodetic altitude) and GC-Alt (geocentric altitude) are my annotations.

[bcoconni]

[...] Geocentric altitude is typically used in orbital mechanics.

I did not know that. Thanks for the information.

[...] whereas geocentric altitude is defined as the height above the ellipsoid surface along a line to the center of the ellipsoid (the radius).

OK, now we have a definition of geocentric altitude for an ellipsoid. Unfortunately JSBSim uses an incorrect definition for the geocentric altitude (the difference between the vehicle radius and the Earth reference radius which is assumed constant).

So something like this? GD-Alt (geodetic altitude) and GC-Alt (geocentric altitude) are my annotations.

The situation is a bit different since the reference point is the vehicle. The contact point P (i.e. the point that intercepts the ellipsoid surface) is different for geocentric and geodetic coordinates.

So even more maths are involved 😢

Right, let's go: the equation of the surface is
$$\left(\frac{x}{a}\right)^2+\left(\frac{y}{b}\right)^2=1$$

The coordinates of the point P are

hence
$$\frac{1}{R^2}=\left(\frac{\cos\psi}{a}\right)^2+\left(\frac{\sin\psi}{b}\right)^2$$

with a bit of trigonometry, we finally get
$$R=\frac{b}{\sqrt{1-e^2\cos^2\psi}}$$

where the eccentricity $e$ is defined by
$$e^2=1-\left(\frac{b}{a}\right)^2$$

and finally the geocentric altitude can be obtained as the difference between the vehicle radius and R.

I'll have to update the code of FGLocation and FGInitialConditions with that.

[bcoconni]

Some further thoughts: according to the definition of the geocentric altitude and knowing that it's only used for orbital mechanics, geocentric altitude only makes sense for a measurement above the sea level, doesn't it ?

On the figure below, I added a green mountain so that it clearly illustrates that AGL might vary greatly depending on the definition that is picked for the geocentric AGL (pink dots below). In addition, I guess AGL is only meaningful for aircraft and irrelevant for orbital mechanics ? So if I am correct, geocentric coordinates are measured above the sea level only, while geodetic coordinates can be measured either above ground level or above sea level.

[seanmcleod]

Just noticed that the latitude values shown in the Wikipedia diagram that I annotated are swapped around. Alpha which is geodetic latitude should be 60 degrees and beta - geocentric latitude should be 66 degrees.

After sketching the annotations I wondered whether we should support (other than for backwards compatibility) the mixing of (lon, lat) coordinates in one datum with an (altitude) in a different datum as you mentioned in one of your examples. In other words if we do support both geodetic and geocentric then maybe require that the altitude component needs to be in the same datum?

In terms of your question regarding geocentric AGL (pink dots) I would imagine the pink dot closer to the vehicle, i.e. on the P-Vehicle dashed line represents the ground level in the geocentric datum. So the geocentric AGL would be the difference between the vehicle's geocentric altitude and the pink dot.

The other pink dot, i.e. the one on the plumb bob line is a ground altitude in the local/vehicle NED datum.

[bcoconni]

Just noticed that the latitude values shown in the Wikipedia diagram that I annotated are swapped around. Alpha which is geodetic latitude should be 60 degrees and beta - geocentric latitude should be 66 degrees.

Hmm... The diagram from Wikipedia seems correct to me. The length CF1 being shorter than AF1, alpha is therefore higher than beta.

After sketching the annotations I wondered whether we should support (other than for backwards compatibility) the mixing of (lon, lat) coordinates in one datum with an (altitude) in a different datum as you mentioned in one of your examples. In other words if we do support both geodetic and geocentric then maybe require that the altitude component needs to be in the same datum?

I tend to agree with you: mixing geocentric and geodetic coordinates should be rejected by JSBSim. Unfortunately we have accumulated some history where the 2 datums have been mixed: this is the case in the example scripts for instance. I have committed some code (16ac672) which allows to specify geocentric latitude with geodetic altitude.

In terms of your question regarding geocentric AGL (pink dots) I would imagine the pink dot closer to the vehicle, i.e. on the P-Vehicle dashed line represents the ground level in the geocentric datum. So the geocentric AGL would be the difference between the vehicle's geocentric altitude and the pink dot.

Drat ! I was expecting an answer like "Sure ! Who cares about geocentric AGL ?" 😄. And unfortunately you picked the dot that requires the more complex calculations but your choice makes sense.

The other pink dot, i.e. the one on the plumb bob line is a ground altitude in the local/vehicle NED datum.

Right, its definition is quite involved but its computation is simpler to implement than the other dot 😄

I am still challenging the relevance of the geocentric AGL however since this altitude is used by orbital mechanics and AFAIK satellites do not require to measure the AGL but rather the ASL.

[seanmcleod]

Hmm... The diagram from Wikipedia seems correct to me.

You're right, I wasn't think straight! 😉

Drat ! I was expecting an answer like "Sure ! Who cares about geocentric AGL ?"

I'm not saying we really have to care about or support the calculation of geocentric AGL, at this point I'm just mentioning what I think the definition of geocentric AGL is.

Now a couple of questions regarding geocentric AGL. Is the aim, if we do decide it's useful, to be able to report geocentric AGL since we support geocentric coordinates in general? So if we support the reporting of geodetic AGL then we should also support the reporting of geocentric AGL?

I haven't looked at the code in detail yet, but it looks like the AGL value for the aircraft is retrieved via the FGGroundCallback interface?

And the FGGroundCallback is passed a FGLocation reference to the vehicle, which internally stores the location in ECEF XYZ cartesian coordinates. The FGDefaultGroundCallback::GetAGLevel() implementation returns a single AGL value, which is the geodetic AGL.

So assuming the following:

If a JSBSim user wanted access to the geocentric AGL without a terrain model, i.e. assuming a spherical earth or ellipsoid, e.g. WGS84 then the ECEF XYZ coordinates would be used to project to a geocentric datum and return the geocentric altitude.

If they wanted access to the geocentric AGL with a terrain model then their FGGroundCallback would need to be able to calculate the terrain elevation above GC-P. If the terrain model data is based on say the WGS84 geodetic datum, e.g. NASA's STRM then they would need to project and interpolate the data in order to calculate a terrain elevation about GC-P. Or I guess simply convert GC-P on the ellipsoid to the equivalent geodetic latitude and then use that to lookup the terrain elevation.

[simfinite]

Glad to see this discussion is alive again. Unfortunately, I am currently not in the position time-wise to contribute much to the solution, but I do want to point out my opinion on GC-AGL vs. GD-AGL:

I think the only time a GC-AGL would be useful in JSBSim (and not considering orbital mechanics), is when assuming a spherical earth. In that case, however, GC-AGL and GD-AGL are identical, right? Thus, the distinction is not required as GC-AGL can be considered a "special case GD-AGL" for a spherical earth.

[bcoconni]

@seanmcleod

Is the aim, if we do decide it's useful, to be able to report geocentric AGL since we support geocentric coordinates in general? So if we support the reporting of geodetic AGL then we should also support the reporting of geocentric AGL?

These are good questions. I already mentioned that I don't think this quantity is relevant but you made a point about the potential inconsistency in JSBSim if we do not support geocentric AGL. I guess that we need a plausible use case to bring that discussion from theoretical highness down to the ground (pun intended 😉) and make a decision. The question is whether there exists a use case for which a missing geocentric AGL computation would make the user's life notably more complicated ?

For example, the case of the geocentric ASL has been settled (as far as I am concerned) when you mentioned orbital mechanics. Indeed the only reference that makes sense for orbital vehicles is the center of the planet since the gravity originates from there.

I haven't looked at the code in detail yet, but it looks like the AGL value for the aircraft is retrieved via the FGGroundCallback interface?

And the FGGroundCallback is passed a FGLocation reference to the vehicle, which internally stores the location in ECEF XYZ cartesian coordinates. The FGDefaultGroundCallback::GetAGLevel() implementation returns a single AGL value, which is the geodetic AGL.

Correct. Terrain is described in terms of geodetic coordinates: the terrain height is the geodetic height of the terrain with respect to the WGS84 ellipsoid measured at the geodetic latitude and longitude passed via the FGLocation instance to FGDefaultGroundCallback::GetAGLevel(). That is how FlightGear proceeds for instance.

If a JSBSim user wanted access to the geocentric AGL without a terrain model, i.e. assuming a spherical earth or ellipsoid, e.g. WGS84 then the ECEF XYZ coordinates would be used to project to a geocentric datum and return the geocentric altitude.

If the terrain definition is missing then the geocentric AGL is trivially equal to the geocentric ASL since the terrain height is zero by default.

If they wanted access to the geocentric AGL with a terrain model then their FGGroundCallback would need to be able to calculate the terrain elevation above GC-P.

Correct. This definition corresponds to the pink dot along the plumb line in my figure above. Its computation is not very expensive because the coordinates of GC-P can be computed quite easily from the equation I derived above. An FGLocation instance can then be constructed from the GC-P coordinates and passed to FGDefaultGroundCallback::GetAGLevel() which would return the AGL of the "plumb line dot".

If the terrain model data is based on say the WGS84 geodetic datum, e.g. NASA's STRM then they would need to project and interpolate the data in order to calculate a terrain elevation about GC-P. Or I guess simply convert GC-P on the ellipsoid to the equivalent geodetic latitude and then use that to lookup the terrain elevation.

The problem with the pink dot that intercepts the terrain along the geocentric radius (see the figure you posted above) is that it is not exactly above GC-P in the geodetic datum. Instead it is slightly offset. The only option (as far as I can see) would be an iterative process that would start from the "plumb line dot" and progress to the seeked intersection. Unfortunately there is no reason the terrain height varies monotonically between the "plumb line dot" and the seeked dot and I guess that will result in some very difficult convergence of the algorithm for a lot of pathological cases.

[bcoconni]

@simfinite

Glad to see this discussion is alive again.

Yes, glad to see it sparks interest even after more than one year idling.

Unfortunately, I am currently not in the position time-wise to contribute much to the solution,

No problem. We all have some real life constraints. You are welcome to bring your point of view to the discussion even if you don't have time to dedicate to coding.

I think the only time a GC-AGL would be useful in JSBSim (and not considering orbital mechanics), is when assuming a spherical earth. In that case, however, GC-AGL and GD-AGL are identical, right? Thus, the distinction is not required as GC-AGL can be considered a "special case GD-AGL" for a spherical earth.

My guts feelings are that your statement is correct. However we are still missing a strong argument to settle the case. Can someone need the geocentric AGL for their specific application ?

[bcoconni]

More food for thoughts: JSBSim is currently using geocentric ASL to compute the atmospheric properties (pressure, temperature, density, humidity rate, etc.). But I think that it should theoretically use geodetic height instead of geocentric ASL because the "J2" gravity varies along the geodetic height. The difference between the 2 altitude definitions is negligible when flying close to the ground but can become measurable beyond stratosphere. Am I over-engineering the case ?

[bcoconni]

FYI, I have asked the Outerra team if they have some usage of the geocentric height and it seems that they are not making any distinction between geodetic and geocentric coordinates. So as far as they are concerned the distinction between geocentric and WGS84 AGL is irrelevant.

I will ask the same question on the FlightGear mailing list and see if we get a different feedback.

[bcoconni]

In the process of updating WGS84 conformance, I have added (commit ecb7952) the ability to modify the planet characteristics with a new XML element <planet>. For example the Earth would be specified with the following WGS84 parameters:

<planet>
  <semimajor_axis unit="M"> 6378137.0 </semimajor_axis>
  <semiminor_axis unit="M"> 6356752.3142 </semiminor_axis>
  <rotation_rate unit="RAD/SEC"> 7292115.8553E-11 </rotation_rate>
  <GM unit="M3/SEC2"> 3986000.9E+8 </GM>
  <J2> 0.00108263 </J2>
</planet>

The idea is to allow JSBSim to simulate flight/orbit/descent on any planet such as the Moon, Mars, etc. This is an idea that has been considered for a long time by the developers as can be seen from the following comments in FGInertial (dating Dec 12, 2008)

jsbsim/src/models/FGInertial.cpp

Lines 67 to 77 in b6e65d4

	// Lunar defaults
	/*
	RotationRate = 0.0000026617;
	GM = 1.7314079E14; // Lunar GM
	RadiusReference = 5702559.05; // Equatorial radius
	C2_0 = 0; // value for the C2,0 coefficient
	J2 = 2.033542482111609E-4; // value for J2
	a = 5702559.05; // semimajor axis length in feet
	b = 5695439.63; // semiminor axis length in feet
	earthPosAngle = 0.0;
	*/

and the existence of the files FGMars.h and FGMars.cpp to model the atmosphere of Mars almost from day one.

[seanmcleod]

I had been meaning to ask how a JSBSim user currently selects between a spherical earth model and the WGS84 ellipsoid model.

Taking a look at the code in FGInertial.cpp it looks to me like the ground callback defaults to using the WGS84 ellipsoid model in the constructor for FGIntertial?

jsbsim/src/models/FGInertial.cpp

Line 73 in ecb7952

GroundCallback.reset(new FGDefaultGroundCallback(a, b));

Then if a planet XML element is specified with different values for the semi-major and semi-minor axes then the ground callback is updated to use these in FGInertial::Load()

jsbsim/src/models/FGInertial.cpp

Line 106 in ecb7952

GroundCallback->SetEllipse(a, b);

So if I'm following the code correctly and haven't missed some other piece it looks like the ground callback defaults to the WGS84 ellipsoid model and if a user wanted to use a spherical earth model they would need to provide a planet XML element with the semi-major and semi-minor axes set to the earth's radius?

[bcoconni]

I had been meaning to ask how a JSBSim user currently selects between a spherical earth model and the WGS84 ellipsoid model.

The short story: it was not possible so far.

Taking a look at the code in FGInertial.cpp it looks to me like the ground callback defaults to using the WGS84 ellipsoid model in the constructor for FGIntertial?

Correct. Even after the commit ecb7952 has been pushed, the default remains the same: Earth is modeled as an oblate spheroid with dimensions from WGS84.

Then if a planet XML element is specified with different values for the semi-major and semi-minor axes then the ground callback is updated to use these in FGInertial::Load()

Yes.

So if I'm following the code correctly and haven't missed some other piece it looks like the ground callback defaults to the WGS84 ellipsoid model and if a user wanted to use a spherical earth model they would need to provide a planet XML element with the semi-major and semi-minor axes set to the earth's radius?

Yes but there is more than the ground callback involved since FGPropagate is also using the semi-major and semi-minor axes:

jsbsim/src/models/FGPropagate.cpp

Lines 122 to 123 in ecb7952

	// For initialization ONLY:
	VState.vLocation.SetEllipse(in.SemiMajor, in.SemiMinor);

and FGInitialConditions also interrogates FGInertialto get the semi-major and semi-minor axes.

jsbsim/src/initialization/FGInitialCondition.cpp

Lines 126 to 129 in ecb7952

	double a = fdmex->GetInertial()->GetSemimajor();
	double b = fdmex->GetInertial()->GetSemiminor();
	position.SetEllipse(a, b);

Previously, FGInertial::GetSemimajor() and FGInertial::GetSemiminor() were still reporting the same values whatever the ground callback was set to. So FGFDMExec::LoadPlanetConstants() has to be called after the planet characteristics have been modified:

jsbsim/src/FGFDMExec.cpp

Lines 721 to 731 in ecb7952

	element = document->FindElement("planet");
	if (element) {
	result = Models[eInertial]->Load(element);
	if (!result) {
	cerr << endl << "Planet element has problems in file " << aircraftCfgFileName << endl;
	return result;
	}
	// Reload the planet constants and re-initialize the models.
	LoadPlanetConstants();
	InitializeModels();
	}

[bcoconni]

More progress with the support of WGS84 in JSBSim: I have recently pushed 2 commits (94fef95 and 679bda6) to address a problem that has been reported a long time ago (SF bug #80 dating from March 17, 2010).

At the moment, the code from these 2 commits is dead code and JSBSim still computes the local frame orientation as it always did. But before proceeding with some further changes, I would like some feedback from the community.

The problem

In JSBSim, the local frame (also known as the NED frame a.k.a. North-East-Down frame - see the definition in the docs) is always oriented so that the down axis points towards the Earth center since the transformation matrix is computed with the geocentric latitude:

jsbsim/src/math/FGLocation.cpp

Lines 301 to 309 in 679bda6

	// Compute the sin/cos values of the latitude
	double sinLat, cosLat;
	if (mRadius == 0.0) {
	sinLat = 0.0;
	cosLat = 1.0;
	} else {
	sinLat = mECLoc(eZ)/mRadius;
	cosLat = rxy/mRadius;
	}

jsbsim/src/math/FGLocation.cpp

Lines 322 to 336 in 679bda6

	// Compute the transform matrices from and to the earth centered frame.
	// See Stevens and Lewis, "Aircraft Control and Simulation", Second Edition,
	// Eqn. 1.4-13, page 40. In Stevens and Lewis notation, this is C_n/e - the
	// orientation of the navigation (local) frame relative to the ECEF frame,
	// and a transformation from ECEF to nav (local) frame.
	mTec2l = { -cosLon*sinLat, -sinLon*sinLat, cosLat,
	-sinLon , cosLon , 0.0 ,
	-cosLon*cosLat, -sinLon*cosLat, -sinLat };
	// In Stevens and Lewis notation, this is C_e/n - the
	// orientation of the ECEF frame relative to the nav (local) frame,
	// and a transformation from nav (local) to ECEF frame.
	mTl2ec = mTec2l.Transposed();

Problem: when the WGS84 gravity is used (which is the default), the down axis as computed by the code above is not parallel to the gravity. This is because the WGS84 gravity does not point towards the Earth center.

In Stevens & Lewis "Aircraft Control and Simulation - 3rd edition", it is specified that the NED frame transformation matrix should be computed with the geodetic latitude so that the down axis be parallel to the spheroid normal (p. 27 for the principles and p. 31 for the formula 1.6-22 of the transformation matrix). They also explain that the gravity computed with the J2 formula (and taking into account the centripetal acceleration) is parallel to the spheroid normal within a few micro-g's (p.33).

This is confirmed by the 2 commits I pushed earlier this week: FGInertial::GetTl2ec() computes the down axis with the J2 gravity and the unit test FGInertialTest::testTransformationMatricesWGS84Earth() checks that the resulting matrix is very close (i.e. within 1e-5 ft/s^2) to the matrix computed with the geodetic latitude.

jsbsim/src/models/FGInertial.cpp

Lines 141 to 167 in 679bda6

	FGMatrix33 FGInertial::GetTl2ec(const FGLocation& location) const
	{
	FGColumnVector3 North, Down, East{-location(eY), location(eX), 0.};
	switch (gravType) {
	case gtStandard:
	{
	Down = location;
	Down *= -1.0;
	}
	break;
	case gtWGS84:
	{
	FGLocation sea_level = location;
	sea_level.SetPositionGeodetic(location.GetLongitude(),
	location.GetGeodLatitudeRad(), 0.0);
	Down = GetGravityJ2(location);
	Down -= vOmegaPlanet*(vOmegaPlanet*sea_level);}
	}
	Down.Normalize();
	East.Normalize();
	North = East*Down;
	return FGMatrix33{North(eX), East(eX), Down(eX),
	North(eY), East(eY), Down(eY),
	North(eZ), 0.0, Down(eZ)};
	}

jsbsim/tests/unit_tests/FGInertialTest.h

Lines 42 to 67 in 679bda6

	void testTransformationMatricesWGS84Earth() {
	FGFDMExec fdmex;
	FGInertial* planet = fdmex.GetInertial();
	fdmex.SetPropertyValue("simulation/gravity-model", 1);
	FGLocation loc;
	FGMatrix33 Tec2l;
	loc.SetEllipse(planet->GetSemimajor(), planet->GetSemiminor());
	for(double lat=-90.; lat <= 90.; lat += 30.) {
	double latitude = lat * degtorad;
	double cosLat = cos(latitude);
	double sinLat = sin(latitude);
	for(double lon=-180.; lon <= 180.; lon += 30.) {
	double longitude = lon * degtorad;
	loc.SetPositionGeodetic(longitude, latitude, 0.0);
	double cosLon = cos(longitude);
	double sinLon = sin(longitude);
	Tec2l = { -cosLon*sinLat, -sinLon*sinLat, cosLat,
	-sinLon , cosLon , 0.0 ,
	-cosLon*cosLat, -sinLon*cosLat, -sinLat };
	TS_ASSERT_MATRIX_EQUALS(planet->GetTec2l(loc), Tec2l);
	TS_ASSERT_MATRIX_EQUALS(planet->GetTl2ec(loc), Tec2l.Transposed());
	}
	}

So what's wrong will you ask me ? The code confirms the statements from Stevens & Lewis and the NED transformation matrices should use the geodetic latitude instead of the geocentric latitude. By the way, that is exactly what has been requested in the SourceForge bug report.

Yes but no

If we just replace the geocentric latitude by the geodetic latitude in FGLocation computations of Tl2ec and Tec2l, we will face 2 other problems:

1. Obviously, the down axis will no longer be aligned with the gravity when the standard gravity will be selected by setting simulation/gravity-model to 0. IMHO the NED frame orientation should depend on the selection of the gravity model (it is the reason why I have coded the new computations of the Tec2l and Tl2ec matrices in the FGInertial class).
2. I am concerned that the fact that the gravity vector being parallel to Earth's spheroid normal is just a fortunate coincidence. If that was always the case then the J2 constant would be a function of the spheroid dimensions only (semi-major and semi-minor lengths) but that is not the case: J2 also depends on the planet mass distribution. So, at least theoretically, there might exist some planets where the gravity vector computed with the J2 formula is not normal to the spheroid. This is a now a problem since commit ecb7952 which allows JSBSim to run simulations for planets other than Earth.

For those 2 reasons, I would like to replace the computation of the NED transformation matrices made by FGLocation (based on the vehicle coordinates) by the computations made by FGInertial (based on the gravity vector).

What do you think ? Any comments or better ideas ?

In other words, the question is should the NED frame be simply tangent to the spheroid ? Or should the definition of the down axis be such that the gravity vector coordinates are always (0, 0, mg) in the NED frame ?

[jonsberndt]

So it seems that the proposed changes would lead to support for both topocentric and topodetic frames? That would be a good thing. Maybe tedious to code, but it would be good. Jon

[seanmcleod]

In other words, the question is should the NED frame be simply tangent to the spheroid ? Or should the definition of the down axis be such that the gravity vector coordinates are always (0, 0, mg) in the NED frame ?

So in picture form from Stevens & Lewis, the difference between the normal to the spheroid versus the normal to the geoid?

Note that the local vertical is defined by the direction in which a plumb-bob hangs and is accurately normal to the geoid. The angle that it makes with the spheroid normal is called the deflection of the vertical and is usually less than 10 arc-s (the largest deflections over the entire Earth are about 1 arc-min).

Now Stevens & Lewis defines NED as:

The local geographic systems have their down axis aligned with the spheroid normal and are oriented North-East-down (ned)

In terms of the definitions for the NED frame in the JSBSim docs the following 3 statements are only all true for a spherical earth model with no geoid model right?

1. plane locally tangent to the Earth’s surface
2. the axis zV points downwards towards the center of the Earth
3. NED convention ensures that the aircraft weight is a force with components (0,0,mg) in the frame

If we go with a NED frame simply tangent to the spheroid what code would be affected?

If a user wanted to compare say some hand-calculated flat-earth calculations (assuming for their hand calculations that they assume the 3 items listed above for the NED frame in the JSBSim documentation currently) then they could achieve that by setting up the following in JSBSim right?

1. Spherical earth via the new planet element
2. Setting simulation/gravity-model = 0

[seanmcleod]

If we go with a NED frame simply tangent to the spheroid what code would be affected?

I don't currently have a good idea of what JSBSim code currently makes use of the NED frame. Came across the following comment while starting to look.

jsbsim/src/models/FGAccelerations.h

Lines 78 to 83 in 2c44e8d

	Special care is taken in the calculations to obtain maximum fidelity in
	JSBSim results. In FGAccelerations, this is obtained by avoiding as much as
	possible the transformations from one frame to another. As a consequence,
	the frames in which the accelerations are primarily evaluated are dictated
	by the frames in which FGPropagate resolves the equations of movement (the
	ECI frame for the translations and the body frame for the rotations).

[bcoconni]

So in picture form from Stevens & Lewis, the difference between the normal to the spheroid versus the normal to the geoid?

Yes, the question is indeed whether we use the local vertical (plumb-bob direction) or the spheroid normal.

In terms of the definitions for the NED frame in the JSBSim docs the following 3 statements are only all true for a spherical earth model with no geoid model right?

1. plane locally tangent to the Earth’s surface
2. the axis zV points downwards towards the center of the Earth
3. NED convention ensures that the aircraft weight is a force with components (0,0,mg) in the frame

Right. These 3 statements can only be simultaneously met for a spherical Earth. The problem with the WGS84 model is that we can meet only one and the question is which one ?

If we go with a NED frame simply tangent to the spheroid what code would be affected?

We would only need to modify the code below so that sinLat and cosLat are the sinus and cosinus of the geodetic latitude (they are currently the sinus and cosinus of the geocentric latitude).

jsbsim/src/math/FGLocation.cpp

Lines 322 to 336 in 679bda6

	// Compute the transform matrices from and to the earth centered frame.
	// See Stevens and Lewis, "Aircraft Control and Simulation", Second Edition,
	// Eqn. 1.4-13, page 40. In Stevens and Lewis notation, this is C_n/e - the
	// orientation of the navigation (local) frame relative to the ECEF frame,
	// and a transformation from ECEF to nav (local) frame.
	mTec2l = { -cosLon*sinLat, -sinLon*sinLat, cosLat,
	-sinLon , cosLon , 0.0 ,
	-cosLon*cosLat, -sinLon*cosLat, -sinLat };
	// In Stevens and Lewis notation, this is C_e/n - the
	// orientation of the ECEF frame relative to the nav (local) frame,
	// and a transformation from nav (local) to ECEF frame.
	mTl2ec = mTec2l.Transposed();

There a number of tests that fail when making that modification however and that needs further investigations.

If a user wanted to compare say some hand-calculated flat-earth calculations (assuming for their hand calculations that they assume the 3 items listed above for the NED frame in the JSBSim documentation currently) then they could achieve that by setting up the following in JSBSim right?

1.Spherical earth via the new planet element
2.Setting simulation/gravity-model = 0

If the user wants to runs simulations around a spherical Earth then yes, the set up would be the one you described above. Alternatively the user could set J2 to 0 in the <planet> element which would have the same effect than setting simulation/gravity-model to 0.

However, flat-Earth cannot be modeled stricto sensu by JSBSim.

[bcoconni]

I don't currently have a good idea of what JSBSim code currently makes use of the NED frame.

There are two main usages:

1. It is used by FGPropagate to compute the attitude angles (pitch, roll, yaw)
2. It is used by FGInitialConditions since initial conditions are mostly specified in the NED frame.

Came across the following comment while starting to look.

Yes I wrote this comment some time ago because the computation of accelerations and propagation had accumulated a number of frame transformations over the years which affected the precision of the calculations. I have submitted a number of patches at the time to reduce the frame transformations to the bare minimum and improve precision.

[bcoconni]

@jonsberndt

So it seems that the proposed changes would lead to support for both topocentric and topodetic frames? That would be a good thing. Maybe tedious to code, but it would be good.

The idea is rather to drop topocentric NED frames and have topodetic NED frames instead. In the event, one would want to use pure topocentric frames, the best option would be to specify a spherical planet with the new <planet> element as @seanmcleod suggested in one of his comments.

[seanmcleod]

It is used by FGPropagate to compute the attitude angles (pitch, roll, yaw)

The angle that it makes with the spheroid normal is called the deflection of the vertical and is usually less than 10 arc-s (the largest deflections over the entire Earth are about 1 arc-min). - Stevens & Lewis

So we should expect to see differences in the pitch angle on that order between a NED frame which has it's down component normal to the spheroid versus one that is normal to the geoid then?

However, flat-Earth cannot be modeled stricto sensu by JSBSim.

Yep. So French, English and Latin at a minimum I see 😉 I had to take Latin as a 3rd language for the 1st two years of high-school.

[bcoconni]

So we should expect to see differences in the pitch angle on that order between a NED frame which has it's down component normal to the spheroid versus one that is normal to the geoid then?

Yes, right.

Following this discussion, the more I think about it, the more I am convinced that I am over-engineering the problem. I guess we will stick with Stevens & Lewis definition of the NED frame (i.e. tangent to the spheroid) and users that mix standard gravity with an oblate planet will have to live with a gravity vector that has a weird direction with respect to the NED frame.

Unless someone has an objection ?

Yep. So French, English and Latin at a minimum I see 😉 I had to take Latin as a 3rd language for the 1st two years of high-school.

Well, that's one of the few Latin expressions I know and since I could not figure out the equivalent English expression I went for that one 😄 Latin expressions are sometimes used in French such as alter ego, I guess it is because the language originates itself from Latin.

[seanmcleod]

users that mix standard gravity with an oblate planet will have to live with a gravity vector that has a weird direction with respect to the NED frame.

Guess we could print a warning if the planet isn't spherical and the user selects standard gravity?

In terms of weird direction in this case are we talking about the difference being the deviation of the normal?

[bcoconni]

Guess we could print a warning if the planet isn't spherical and the user selects standard gravity?

I don't know. We already have an issue #140 opened about having too much warnings.

In terms of weird direction in this case are we talking about the difference being the deviation of the normal?

That's correct. 11.5 arc-min is not much so how low the deviation D should be before issuing the warning ? Do you have some idea about what the criteria could be ?

[seanmcleod]

I wasn't thinking of having the warning printed based on a particular deviation value, rather a general warning that the user is using a non-standard setup in terms of standard gravity combined with a non-spherical planet.

I was also thinking in terms of a warning message on startup. If you were to look at the deviation of the normal and only display a message based on some cut-off value then there may be no warning on startup but only later on as the aircraft flies north/south right?

In terms of the too many warnings I guess we could look at having some warning levels that the user can suppress, or disable specific warnings?

[bcoconni]

I wasn't thinking of having the warning printed based on a particular deviation value, rather a general warning that the user is using a non-standard setup in terms of standard gravity combined with a non-spherical planet.

Well, what would be the criterion to decide if the planet is non-spherical ? As soon as the semi-major axis a is different than the semi-minor axis b not matter how much different they are ? Or if (a-b)/a is below a threshold ?

If you were to look at the deviation of the normal and only display a message based on some cut-off value then there may be no warning on startup but only later on as the aircraft flies north/south right?

Since we know that the maximum deviation occurs at a latitude of 45 degrees, this can be evaluated as soon as simulation/gravity-model is set to 0. But the question remains the same than above: what would be the threshold above which the warning should be issued ?

My personal preference would be the 2nd option. What do you think ?

In terms of the too many warnings I guess we could look at having some warning levels that the user can suppress, or disable specific warnings?

Yes, I think we should migrate to FlightGear/SimGear's log system SG_LOG which allows that kind of settings. Yet another item on our TODO list 😉

[seanmcleod]

My thinking is not to try and determine a threshold related to accuracy and to warn the user then. Rather only to warn the user for some specific combinations that may indicate they've made a mistake in the combination they've chosen, i.e.

Standard gravity and non-spherical planet (a != b) - thought they were using a spherical model and standard gravity.

EGM96 J2 gravity and spherical planet (a == b) - thought they were using WGS84 spheroid model EGM96 J2 gravity.

If they did mean to use the specific combination then some option to suppress the warning.

In taking a look at some of the recent code for the planet element etc. it looks like a user can end up with the specific EGM96 J2 value for earth if they forget to supply it as part of the planet element?

jsbsim/src/models/FGInertial.cpp

Line 58 in ecb7952

J2 = 1.08262982E-03; // WGS84 value for J2

In other words should we display an error if a planet is specified with axes that don't match WGS84 and no J2 value is supplied and the J2 gravity model is selected?

[bcoconni]

@seanmcleod
I have submitted the PR #318 with a proposal implementation of the warnings messages that you suggested (plus some other to cover all the combinations 😄). Please review and let me know what you think.
Also since you're a native English speaker, by all means do not hesitate to rephrase the warning messages !

[bcoconni]

For the record, I just pushed the commit 70a327f which converts the waypoint calculations to WGS84.

FGLocation has 2 methods GetDistanceTo and GetHeadingTo which allow to compute respectively the geodesic distance s12 and the azimuth α1 (i.e. the heading) between 2 points A and B located on the WGS84 ellipsoid.

Previously, these 2 methods were using the Haversine formulas which are only valid for a spherical planet. Unfortunately the computation of the geodesic distance and azimuth on an ellipsoid is quite complex so, rather than developing our own code from scratch, I have used Charles F. F. Karney's library GeographicLib. GeographicLib is licensed under the MIT/X11 license so I have included the source code from GeographicLib that is relevant to the computation of geodesic distances and headings. That way, JSBSim still has no external dependencies.

The author of GeographicLib has conducted quite thorough investigations on the topic as can be seen from the very impressive bibliography that he has collected. He is also maintaining his library for several years now so I am quite confident that we can rely on his library and save ourselves the trouble of maintaining such a piece of complex and very specialized software.

This also updates the <waypoint> flight control in the process since it relies on FGLocation to process the distance and heading between two points.

In case some users have the GeographicLib library already installed on their computer, I plan to update our CMake build process to be able to link with an external GeographicLib library when there is one, otherwise our copy of the GeographicLib source code will be compiled.

As far as I can tell this was the last bit of JSBSim that was not conformant with WGS84 so, unless someone objects, I plan to close this issue which is now more than 2 years old.

[simfinite]

Thanks a lot for working on this. I think the WGS84 support will really push JSBSim's applicability for real-world scenarios where aeronautical navigation is required.

[bcoconni]

@simfinite You're welcome !

So I am closing this issue. Would any of you find an error, a bug or a feature that still needs to be converted to WGS84, please open a new issue.