Add FGNullAtmosphere and extend LoadPlanet atmosphere factory (JSBSim-Team#1462)

Co-authored-by: Bertrand Coconnier <bcoconni@users.noreply.github.com>

Author: AlexanderRex

src/FGFDMExec.cpp

@@ -47,6 +47,7 @@ INCLUDES
 #include "FGFDMExec.h"
 #include "models/atmosphere/FGStandardAtmosphere.h"
 #include "models/atmosphere/FGMSIS.h"
+#include "models/atmosphere/FGMars.h"
 #include "models/atmosphere/FGWinds.h"
 #include "models/FGFCS.h"
 #include "models/FGPropulsion.h"
@@ -63,6 +64,7 @@ INCLUDES
 #include "initialization/FGLinearization.h"
 #include "input_output/FGScript.h"
 #include "input_output/FGXMLFileRead.h"
+#include "input_output/string_utilities.h"
 #include "initialization/FGInitialCondition.h"
 #include "input_output/FGLog.h"
 
@@ -788,10 +790,20 @@ bool FGFDMExec::LoadPlanet(Element* element)
     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
+      to_lower(model);
+
+      if (model == "msis" || model == "mars") {
         instance->Unbind(Models[eAtmosphere]);
-        Models[eAtmosphere] = std::make_shared<FGMSIS>(this);
+      }
+
+      std::shared_ptr<FGAtmosphere> newAtmosphere;
+      if (model == "msis")
+        newAtmosphere = std::make_shared<FGMSIS>(this);
+      else if (model == "mars")
+        newAtmosphere = std::make_shared<FGMars>(this);
+
+      if (newAtmosphere) {
+        Models[eAtmosphere] = newAtmosphere;
         Atmosphere = static_cast<FGAtmosphere*>(Models[eAtmosphere].get());
 
         // Model initialization sequence

src/models/atmosphere/CMakeLists.txt

@@ -1,10 +1,12 @@
 set(SOURCES FGMSIS.cpp
+            FGMars.cpp
             FGStandardAtmosphere.cpp
             FGWinds.cpp
             MSIS/nrlmsise-00.c
             MSIS/nrlmsise-00_data.c)
 
 set(HEADERS FGMSIS.h
+            FGMars.h
             FGStandardAtmosphere.h
             FGWinds.h
             MSIS/nrlmsise-00.h)

src/models/atmosphere/FGMars.cpp

@@ -61,7 +61,6 @@ FGMars::FGMars(FGFDMExec* fdmex) : FGAtmosphere(fdmex)
   Name = "FGMars";
   Reng = 53.5 * 44.01;
 
-  bind();
   Debug(0);
 }
 
@@ -74,9 +73,9 @@ void FGMars::Calculate(double altitude)
   // LIMIT the temperatures so they do not descend below absolute zero.
 
   if (altitude < 22960.0) {
-    Temperature = -25.68 - 0.000548*altitude; // Deg Fahrenheit
+    Temperature = -25.68 - 0.000548*altitude + 459.67; // Convert F to Rankine
   } else {
-    Temperature = -10.34 - 0.001217*altitude; // Deg Fahrenheit
+    Temperature = -10.34 - 0.001217*altitude + 459.67; // Convert F to Rankine
   }
   Pressure = 14.62*exp(-0.00003*altitude); // psf - 14.62 psf =~ 7 millibars
   Density = Pressure/(Reng*Temperature); // slugs/ft^3 (needs deg R. as input

src/models/atmosphere/FGMars.h

@@ -63,10 +63,13 @@ class FGMars : public FGAtmosphere {
   /// Constructor
   FGMars(FGFDMExec*);
 
-private:
-  void Calculate(double altitude);
+  double GetTemperature(double altitude) const override { return Temperature; }
+  double GetPressure(double altitude) const override { return Pressure; }
+  void SetTemperature(double t, double h, eTemperature unit) override {}
 
-  void Debug(int from);
+private:
+  void Calculate(double altitude) override;
+  void Debug(int from) override;
 };
 
 } // namespace JSBSim

tests/TestPlanet.py

@@ -75,6 +75,49 @@ def test_load_MSIS_atmosphere(self):
         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_mars_atmosphere(self):
+        # Mars atmosphere via <planet><atmosphere model="Mars"/></planet>.
+        # Reference values are the closed-form output of FGMars::Calculate at
+        # altitude = 0:
+        #   T = -25.68 + 459.67 = 433.99 R         (~241.1 K)
+        #   P = 14.62 psf                          (~7 mbar)
+        #   rho = P / (Reng * T), Reng = 53.5*44.01 (CO2)
+        tripod = FlightModel(self, 'tripod')
+        mars_file = self.sandbox.path_to_jsbsim_file('tests/mars.xml')
+        tripod.include_planet_test_file(mars_file)
+        self.fdm = tripod.start()
+        self.fdm['ic/h-agl-ft'] = 0.0
+        self.fdm['ic/long-gc-deg'] = 0.0
+        self.fdm['ic/lat-geod-deg'] = 0.0
+        self.fdm.run_ic()
+
+        # Mars equatorial radius = 3396.2 km
+        self.assertAlmostEqual(self.fdm['metrics/terrain-radius']*0.3048/3396200, 1.0)
+        # Surface gravity ~3.71-3.72 m/s^2 (JSBSim uses oblate-spheroid gravity,
+        # not point-mass GM/R^2, so the value is a few mGal above textbook).
+        self.assertAlmostEqual(self.fdm['accelerations/gravity-ft_sec2']*0.3048, 3.72, delta=5E-2)
+
+        self.assertAlmostEqual(self.fdm['atmosphere/T-R'], 433.99, delta=1E-2)
+        self.assertAlmostEqual(self.fdm['atmosphere/P-psf'], 14.62, delta=1E-2)
+        # rho = 14.62 / (53.5*44.01 * 433.99) ~ 1.4308e-5 slugs/ft^3
+        self.assertAlmostEqual(self.fdm['atmosphere/rho-slugs_ft3'], 1.4308e-5, delta=1E-8)
+
+    def test_load_Mars_atmosphere(self):
+        # Same as above but using FGFDMExec::LoadPlanet at runtime instead of
+        # including the planet file at FDM construction time.
+        tripod = FlightModel(self, 'tripod')
+        mars_file = self.sandbox.path_to_jsbsim_file('tests/mars.xml')
+        self.fdm = tripod.start()
+        self.fdm.load_planet(mars_file, False)
+        self.fdm['ic/h-agl-ft'] = 0.0
+        self.fdm['ic/long-gc-deg'] = 0.0
+        self.fdm['ic/lat-geod-deg'] = 0.0
+        self.fdm.run_ic()
+
+        self.assertAlmostEqual(self.fdm['atmosphere/T-R'], 433.99, delta=1E-2)
+        self.assertAlmostEqual(self.fdm['atmosphere/P-psf'], 14.62, delta=1E-2)
+        self.assertAlmostEqual(self.fdm['atmosphere/rho-slugs_ft3'], 1.4308e-5, delta=1E-8)
+
     def test_planet_geographic_error1(self):
         # Check that a negative equatorial radius raises an exception
         tripod = FlightModel(self, 'tripod')

tests/mars.xml

@@ -0,0 +1,8 @@
+<planet name="Mars">
+  <equatorial_radius unit="KM">3396.2</equatorial_radius>
+  <polar_radius unit="KM">3376.2</polar_radius>
+  <rotation_rate unit="RAD/SEC">7.0882E-5</rotation_rate>
+  <GM unit="M3/SEC2">4.2828372E13</GM>
+  <J2>1.9605E-3</J2>
+  <atmosphere model="Mars"/>
+</planet>