Enable the usage of the NRLMSIS00 atmospheric model. (#941)

src/FGFDMExec.cpp

@@ -46,6 +46,7 @@ INCLUDES
 
 #include "FGFDMExec.h"
 #include "models/atmosphere/FGStandardAtmosphere.h"
+#include "models/atmosphere/FGMSIS.h"
 #include "models/atmosphere/FGWinds.h"
 #include "models/FGFCS.h"
 #include "models/FGPropulsion.h"
@@ -451,7 +452,9 @@ void FGFDMExec::LoadInputs(unsigned int idx)
     Inertial->in.Position      = Propagate->GetLocation();
     break;
   case eAtmosphere:
-    Atmosphere->in.altitudeASL = Propagate->GetAltitudeASL();
+    Atmosphere->in.altitudeASL     = Propagate->GetAltitudeASL();
+    Atmosphere->in.GeodLatitudeDeg = Propagate->GetGeodLatitudeDeg();
+    Atmosphere->in.LongitudeDeg    = Propagate->GetLongitudeDeg();
     break;
   case eWinds:
     Winds->in.AltitudeASL      = Propagate->GetAltitudeASL();
@@ -762,7 +765,7 @@ bool FGFDMExec::LoadPlanet(const SGPath& PlanetPath, bool useAircraftPath)
 
   if (document->GetName() != "planet") {
     stringstream s;
-    s << "File: " << PlanetFileName << " is not a reset file.";
+    s << "File: " << PlanetFileName << " is not a planet file.";
     cerr << s.str() << endl;
     throw BaseException(s.str());
   }
@@ -786,6 +789,28 @@ bool FGFDMExec::LoadPlanet(Element* element)
     LoadPlanetConstants();
     IC->InitializeIC();
     InitializeModels();
+
+    // Process the atmosphere element. This element is OPTIONAL.
+    Element* atm_element = element->FindElement("atmosphere");
+    if (atm_element && atm_element->HasAttribute("model")) {
+      string model = atm_element->GetAttributeValue("model");
+      if (model == "MSIS") {
+        // Replace the existing atmosphere model
+        instance->Unbind(Models[eAtmosphere].get());
+        Models[eAtmosphere] = std::make_shared<FGMSIS>(this);
+        Atmosphere = static_cast<FGAtmosphere*>(Models[eAtmosphere].get());
+
+        // Model initialization sequence
+        LoadInputs(eAtmosphere);
+        Atmosphere->InitModel();
+        result = Atmosphere->Load(atm_element);
+        if (!result) {
+          cerr << endl << "Incorrect definition of <atmosphere>." << endl;
+          return result;
+        }
+        InitializeModels();
+      }
+    }
   }
 
   return result;

tests/MSIS.xml

@@ -0,0 +1,6 @@
+<planet name="Earth">
+  <atmosphere model="MSIS">
+    <day>172</day>    <!--  Jun 21st  -->
+    <utc>29000</utc>  <!-- 08:03:20AM -->
+  </atmosphere>
+</planet>

tests/TestPlanet.py

@@ -61,6 +61,20 @@ def test_load_planet(self):
 
         self.assertAlmostEqual(self.fdm['metrics/terrain-radius']*0.3048/1736000, 1.0)
 
+    def test_load_MSIS_atmosphere(self):
+        tripod = FlightModel(self, 'tripod')
+        MSIS_file = self.sandbox.path_to_jsbsim_file('tests/MSIS.xml')
+        self.fdm = tripod.start()
+        self.fdm.load_planet(MSIS_file, False)
+        self.fdm['ic/h-sl-ft'] = 0.0
+        self.fdm['ic/long-gc-deg'] = -70.0
+        self.fdm['ic/lat-geod-deg'] = 60.0
+        self.fdm.run_ic()
+
+        self.assertAlmostEqual(self.fdm['atmosphere/T-R']*5/9, 281.46476, delta=1E-5)
+        self.assertAlmostEqual(self.fdm['atmosphere/rho-slugs_ft3']/0.001940318, 1.263428, delta=1E-6)
+        self.assertAlmostEqual(self.fdm['atmosphere/P-psf'], 2132.294, delta=1E-3)
+
     def test_planet_geographic_error1(self):
         # Check that a negative equatorial radius raises an exception
         tripod = FlightModel(self, 'tripod')