Improvements to the density altitude computation about jsbsim HOT 12 CLOSED

@djlinse thanks for the comments. Yep the equation which I grabbed from a much earlier version of JSBSim is only valid up to 36089ft since it uses a single lapse rate. It seems to be the common/standard formula you come across on Wikipedia etc. for density altitude up to 36089ft assuming dry air.

But until we implement one that is valid at higher altitudes we should document it's current limitation. I'll take a look at the link you posted.

I'll also look into the code structure suggestions.

from jsbsim.

@bcoconni @djlinse on some further studying of the atmosphere code I noticed a potential issue with a difference between the geometric altitude (GeoMet) and the geopotential altitude (GeoPot).

Take a look at the following code from the test snippet that @bcoconni submitted.

jsbSim_density_alt = fdm['atmosphere/density-altitude']
# Check that when setting the aircraft at the altitude jsbSim_density_alt with delta-T == 0.0
# we get the same density than at the pressure altitude with delta-T != 0.0
density_ref = fdm['atmosphere/rho-slugs_ft3']

fdm['atmosphere/delta-T'] = 0.0
fdm['ic/h-agl-ft'] = jsbSim_density_alt
fdm.run_ic()

self.assertAlmostEqual(density_ref, fdm['atmosphere/rho-slugs_ft3'], delta=1e-5)

The altitude passed to fdm['ic/h-agl-ft'] is assumed to be a GeoMet altitude and so the atmosphere code converts it to a GeoPot altitude before making use of the atmosphere tables etc.

And in the current implementation of the DensityAltitude() function we return a GeoPot altitude, which the test then passes back in via fdm['ic/h-agl-ft'] = jsbSim_density_alt as a GeoMet altitude.

If I return a GeoMet altitude from the DensityAltitude() function then the error differences improve by 1-2 orders of magnitude from 10,000ft upwards, no real difference for the sea-level test.

from jsbsim.

@bcoconni @djlinse I've submitted a new pull request 56 to address these issues except for a density altitude calculation above the Troposphere, will look into that next.

from jsbsim.

Sean,

In general, I like it. Not surprising, I suppose, given you've changed what I asked for! A few minor, technical suggestions:

1. "Troposphere", not "toposphere", in a number of function and variable names and comments.
2. A copy/paste left you with double semicolons (;;) on some of the calls to GeopotentialAltitude().
3. For this bit:

// Earth radius in ft as defined for ISA 1976
const double FGStandardAtmosphere::EarthRadius = 20855531.5;

I might say something like this (taking from 1976 ISA)

// Effective radius of the earth at a specific latitude per ISA 1976 (converted to ft)
// r0 = 6356766 m
const double FGStandardAtmosphere::EarthRadius = 6356766.0/0.3048;

4. The comment before DAToposphereExponent.

Per Eqn 12 in https://wahiduddin.net/calc/density_altitude.htm, this exponent is LR/(g0M - L*R). It is unitless, in the end, but does have the lapse rate, L, in it, so there are a couple of choices.

a. Easiest, just accept the ISA standard values where L is defined in metric units as -6.5 K/km (for the troposphere), and all of the other values are also defined values in the ISA standard. So ~0.235... is calculated with those values one time.

b. If you think it is important (I personally don't), you could recalculate the exponent using the current (possibly non-standard) lapse rate. Just convert that lapse rate to metric and use the ISA values for R, g0, and M.

Actually, b. is wrong, in that we're trying to figure out the altitude in the standard atmosphere. So, either in comment or calculation just say how the value is determined.

from jsbsim.

Dennis thanks for the feedback. I've submitted a new commit to address your feedback.

The value for the DATroposphereExponent had been calculated offline using the SI units for R, g0 and M and the ISA 1976 lapse rate L for the Troposphere. I didn't want to introduce SI units for g0, R, M etc. just to compute this exponent, so I was looking to use the existing British units in FGJSBBase, or to leave it as a pre-computed value from SI units etc.

from jsbsim.

Sean, it looks fine to me. If it were a particularly important constant, I might suggest just having a longer comment that includes the SI values and says that they come directly from the ISA document. This is not really such a constant.

from jsbsim.

We must also check that the density altitude you get from there is consistent with the density computations made in the C++ class FGStandardAtmosphere

If the test TestDensityAltitude.py is modified as below

        for pressure_alt, delta_T, density_alt in reference_data:

            fdm['ic/h-sl-ft'] = pressure_alt
            fdm['atmosphere/delta-T'] = delta_T
            fdm.run_ic()

            jsbSim_density_alt = fdm['atmosphere/density-altitude']

            self.assertAlmostEqual(density_alt, jsbSim_density_alt, delta=1e-7)

            # Check that when setting the aircraft at the altitude jsbSim_density_alt with delta-T == 0.0
            # we get the same density than at the pressure altitude with delta-T != 0.0
            density_ref = fdm['atmosphere/rho-slugs_ft3']

            fdm['ic/h-sl-ft'] = jsbSim_density_alt
            fdm['atmosphere/delta-T'] = 0.0
            fdm.run_ic()

            self.assertAlmostEqual(density_ref, fdm['atmosphere/rho-slugs_ft3'])

then the test fails with the following error

 ======================================================================
 FAIL: test_densityaltitude (__main__.TestDensityAltitude)
 ----------------------------------------------------------------------
 Traceback (most recent call last):
   File "TestDensityAltitude.py", line 61, in test_densityaltitude
     self.assertAlmostEqual(density_ref, fdm['atmosphere/rho-slugs_ft3'])
 AssertionError: 0.002507425951046395 != 0.0025391558128942003 within 7 place

which means that the formula used to compute the density altitude has an inconsistency with the other formulas in FGStandardAtmosphere that computes the density.

We must find why.

Note: I think we should use the fdm['ic/h-sl-ft'] (aka altitude above sea level) in the TestDensityAltitude rather than fdm['ic/h-agl-ft'] (aka altitude above ground level) since the output of FGStandardAtmosphere::GetDensityAltitude() is an ASL. I have modified the code snippet above accordingly.

from jsbsim.

@bcoconni I'll update the test, just making the requested changes to the C++ code first.

In terms of the consistency with C++ density calculations take a look at my earlier post above in this thread where the difference improved to be on the order of 1e-8 to 1-e9 for all the tests except for the sea level -27R which is 3.17e-5.

from jsbsim.

In terms of the consistency with C++ density calculations take a look at my earlier post above in this thread where the difference improved to be on the order of 1e-8 to 1-e9 for all the tests except for the sea level -27R which is 3.17e-5.

Yes, sorry. I was focused on copying here the comment I made on your PR #38 for the recored. But did not make an update with regard to your progress, so this breaks the thread of discussion. Sorry about that.

I noticed the improvement you got now that the geopotential altitude is converted to geometric altitude and that's great 👍 . We are almost there.

However there remains this difference at sea level for delta-T == -27°F that we cannot explain. I am trying to investigate on my side as well. With an additional pair of eyes we should be able to squash this bug 😃

from jsbsim.

@bcoconni no problem. I've spent a lot more time refreshing my knowledge of the ISA than I expected after initially copying in the density altitude calculation from a much older commit of JSBSim, but it has been educational ;-)

In the meantime I've also got some code ready for a separate issue in terms of calculating the pressure altitude for a non-standard ISA day.

from jsbsim.

Caught it ! The code in FGStandardAtmosphere is not designed for negative altitudes (which is the output for GetDensityAltitude(0.0) for any negative delta-T

🐛 💥 🔫

from jsbsim.

OK the bug is fixed (commit f0e7d64)

from jsbsim.