Move the Mach/VCAS conversion methods to FGAtmosphere (#881)

* Remove `static` methods and members
* Replace magic constants with their expression
* Use altitude to compute VCAS or Mach instead of ambient pressure
* Improved unit test

src/FGFDMExec.cpp

@@ -470,6 +470,7 @@ void FGFDMExec::LoadInputs(unsigned int idx)
  Auxiliary->in.SoundSpeed = Atmosphere->GetSoundSpeed();
  Auxiliary->in.KinematicViscosity = Atmosphere->GetKinematicViscosity();
  Auxiliary->in.DistanceAGL = Propagate->GetDistanceAGL();
+ Auxiliary->in.AltitudeASL = Propagate->GetAltitudeASL();
  Auxiliary->in.Mass = MassBalance->GetMass();
  Auxiliary->in.Tl2b = Propagate->GetTl2b();
  Auxiliary->in.Tb2l = Propagate->GetTb2l();

src/FGJSBBase.cpp

@@ -102,72 +102,6 @@ string FGJSBBase::CreateIndexedPropertyName(const string& Property, int index)
  buf << Property << '[' << index << ']';
  return buf.str();
 }
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-double FGJSBBase::PitotTotalPressure(double mach, double p)
-{
- if (mach < 0) return p;
- if (mach < 1) //calculate total pressure assuming isentropic flow
- return p*pow((1 + 0.2*mach*mach),3.5);
- else {
- // shock in front of pitot tube, we'll assume its normal and use
- // the Rayleigh Pitot Tube Formula, i.e. the ratio of total
- // pressure behind the shock to the static pressure in front of
- // the normal shock assumption should not be a bad one -- most supersonic
- // aircraft place the pitot probe out front so that it is the forward
- // most point on the aircraft. The real shock would, of course, take
- // on something like the shape of a rounded-off cone but, here again,
- // the assumption should be good since the opening of the pitot probe
- // is very small and, therefore, the effects of the shock curvature
- // should be small as well. AFAIK, this approach is fairly well accepted
- // within the aerospace community
-
- // The denominator below is zero for Mach ~ 0.38, for which
- // we'll never be here, so we're safe
-
- return p*166.92158009316827*pow(mach,7.0)/pow(7*mach*mach-1,2.5);
- }
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-// Based on the formulas in the US Air Force Aircraft Performance Flight Testing 
-// Manual (AFFTC-TIH-99-01). In particular sections 4.6 to 4.8.
-
-double FGJSBBase::MachFromImpactPressure(double qc, double p)
-{
- double A = qc / p + 1;
- double M = sqrt(5.0*(pow(A, 1. / 3.5) - 1)); // Equation (4.12)
-
- if (M > 1.0)
- for (unsigned int i = 0; i<10; i++)
- M = 0.8812848543473311*sqrt(A*pow(1 - 1.0 / (7.0*M*M), 2.5)); // Equation (4.17)
-
- return M;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-double FGJSBBase::VcalibratedFromMach(double mach, double p)
-{
- double asl = FGAtmosphere::StdDaySLsoundspeed;
- double psl = FGAtmosphere::StdDaySLpressure;
- double qc = PitotTotalPressure(mach, p) - p;
-
- return asl * MachFromImpactPressure(qc, psl); 
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-double FGJSBBase::MachFromVcalibrated(double vcas, double p)
-{
- double asl = FGAtmosphere::StdDaySLsoundspeed;
- double psl = FGAtmosphere::StdDaySLpressure;
- double qc = PitotTotalPressure(vcas / asl, psl) - psl;
-
- return MachFromImpactPressure(qc, p);
-}
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 } // namespace JSBSim

src/FGJSBBase.h

@@ -250,41 +250,6 @@ class JSBSIM_API FGJSBBase {
  return measure*0.3048;
  }
 
- /** Compute the total pressure in front of the Pitot tube. It uses the
- * Rayleigh formula for supersonic speeds (See "Introduction to Aerodynamics
- * of a Compressible Fluid - H.W. Liepmann, A.E. Puckett - Wiley & sons
- * (1947)" §5.4 pp 75-80)
- * @param mach The Mach number
- * @param p Pressure in psf
- * @return The total pressure in front of the Pitot tube in psf */
- static double PitotTotalPressure(double mach, double p);
-
- /** Compute the Mach number from the differential pressure (qc) and the
- * static pressure. Based on the formulas in the US Air Force Aircraft
- * Performance Flight Testing Manual (AFFTC-TIH-99-01).
- * @param qc The differential/impact pressure
- * @param p Pressure in psf
- * @return The Mach number */
- static double MachFromImpactPressure(double qc, double p);
-
- /** Calculate the calibrated airspeed from the Mach number. Based on the
- * formulas in the US Air Force Aircraft Performance Flight Testing
- * Manual (AFFTC-TIH-99-01).
- * @param mach The Mach number
- * @param p Pressure in psf
- * @return The calibrated airspeed (CAS) in ft/s
- * */
- static double VcalibratedFromMach(double mach, double p);
-
- /** Calculate the Mach number from the calibrated airspeed.Based on the
- * formulas in the US Air Force Aircraft Performance Flight Testing
- * Manual (AFFTC-TIH-99-01).
- * @param vcas The calibrated airspeed (CAS) in ft/s
- * @param p Pressure in psf
- * @return The Mach number
- * */
- static double MachFromVcalibrated(double vcas, double p);
-
  /** Finite precision comparison.
  @param a first value to compare
  @param b second value to compare

src/initialization/FGInitialCondition.cpp

@@ -174,8 +174,7 @@ void FGInitialCondition::SetMachIC(double mach)
 void FGInitialCondition::SetVcalibratedKtsIC(double vcas)
 {
  double altitudeASL = GetAltitudeASLFtIC();
- double pressure = Atmosphere->GetPressure(altitudeASL);
- double mach = MachFromVcalibrated(fabs(vcas)*ktstofps, pressure);
+ double mach = Atmosphere->MachFromVcalibrated(fabs(vcas)*ktstofps, altitudeASL);
  double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
 
  SetVtrueFpsIC(mach * soundSpeed);
@@ -677,13 +676,12 @@ double FGInitialCondition::GetTerrainElevationFtIC(void) const
 void FGInitialCondition::SetAltitudeAGLFtIC(double agl)
 {
  double altitudeASL = GetAltitudeASLFtIC();
- double pressure = Atmosphere->GetPressure(altitudeASL);
  double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
  double rho = Atmosphere->GetDensity(altitudeASL);
  double rhoSL = Atmosphere->GetDensitySL();
 
  double mach0 = vt / soundSpeed;
- double vc0 = VcalibratedFromMach(mach0, pressure);
+ double vc0 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL);
  double ve0 = vt * sqrt(rho/rhoSL);
 
  switch(lastLatitudeSet) {
@@ -721,11 +719,10 @@ void FGInitialCondition::SetAltitudeAGLFtIC(double agl)
  altitudeASL = GetAltitudeASLFtIC();
  soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
  rho = Atmosphere->GetDensity(altitudeASL);
- pressure = Atmosphere->GetPressure(altitudeASL);
 
  switch(lastSpeedSet) {
  case setvc:
- mach0 = MachFromVcalibrated(vc0, pressure);
+ mach0 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL);
  SetVtrueFpsIC(mach0 * soundSpeed);
  break;
  case setmach:
@@ -749,13 +746,12 @@ void FGInitialCondition::SetAltitudeAGLFtIC(double agl)
 void FGInitialCondition::SetAltitudeASLFtIC(double alt)
 {
  double altitudeASL = GetAltitudeASLFtIC();
- double pressure = Atmosphere->GetPressure(altitudeASL);
  double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
  double rho = Atmosphere->GetDensity(altitudeASL);
  double rhoSL = Atmosphere->GetDensitySL();
 
  double mach0 = vt / soundSpeed;
- double vc0 = VcalibratedFromMach(mach0, pressure);
+ double vc0 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL);
  double ve0 = vt * sqrt(rho/rhoSL);
 
  switch(lastLatitudeSet) {
@@ -825,11 +821,10 @@ void FGInitialCondition::SetAltitudeASLFtIC(double alt)
  altitudeASL = position.GetGeodAltitude();
  soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
  rho = Atmosphere->GetDensity(altitudeASL);
- pressure = Atmosphere->GetPressure(altitudeASL);
 
  switch(lastSpeedSet) {
  case setvc:
- mach0 = MachFromVcalibrated(vc0, pressure);
+ mach0 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL);
  SetVtrueFpsIC(mach0 * soundSpeed);
  break;
  case setmach:
@@ -959,11 +954,10 @@ double FGInitialCondition::GetBodyWindFpsIC(int idx) const
 double FGInitialCondition::GetVcalibratedKtsIC(void) const
 {
  double altitudeASL = GetAltitudeASLFtIC();
- double pressure = Atmosphere->GetPressure(altitudeASL);
  double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
  double mach = vt / soundSpeed;
 
- return fpstokts * VcalibratedFromMach(mach, pressure);
+ return fpstokts * Atmosphere->VcalibratedFromMach(mach, altitudeASL);
 }
 
 //******************************************************************************

src/models/FGAtmosphere.cpp

@@ -51,11 +51,9 @@ namespace JSBSim {
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-// Atmosphere constants in British units converted from the SI values specified in the
-// ISA document - https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf
-const double FGAtmosphere::StdDaySLsoundspeed = sqrt(SHRatio*Reng0*StdDaySLtemperature);
-
-FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
+FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex)
+ : FGModel(fdmex),
+ StdDaySLsoundspeed(sqrt(SHRatio*Reng0*StdDaySLtemperature))
 {
  Name = "FGAtmosphere";
 
@@ -283,6 +281,8 @@ double FGAtmosphere::ConvertToPSF(double p, ePressure unit) const
  return targetPressure;
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 double FGAtmosphere::ConvertFromPSF(double p, ePressure unit) const
 {
  double targetPressure=0; // Pressure
@@ -309,6 +309,79 @@ double FGAtmosphere::ConvertFromPSF(double p, ePressure unit) const
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+double FGAtmosphere::PitotTotalPressure(double mach, double p) const
+{
+ constexpr double a = (SHRatio-1.0) / 2.0;
+ constexpr double b = SHRatio / (SHRatio-1.0);
+ constexpr double c = 2.0*b;
+ constexpr double d = 1.0 / (SHRatio-1.0);
+ const double coeff = pow(0.5*(SHRatio+1.0), b)*pow((SHRatio+1.0)/(SHRatio-1.0), d);
+
+ if (mach < 0) return p;
+ if (mach < 1) //calculate total pressure assuming isentropic flow
+ return p*pow((1.0 + a*mach*mach), b);
+ else {
+ // Shock in front of pitot tube, we'll assume its normal and use the
+ // Rayleigh Pitot Tube Formula, i.e. the ratio of total pressure behind the
+ // shock to the static pressure in front of the normal shock assumption
+ // should not be a bad one -- most supersonic aircraft place the pitot probe
+ // out front so that it is the forward most point on the aircraft.
+ // The real shock would, of course, take on something like the shape of a
+ // rounded-off cone but, here again, the assumption should be good since the
+ // opening of the pitot probe is very small and, therefore, the effects of
+ // the shock curvature should be small as well. AFAIK, this approach is
+ // fairly well accepted within the aerospace community
+
+ // The denominator below is zero for Mach ~ 0.38, for which
+ // we'll never be here, so we're safe
+
+ return p*coeff*pow(mach, c)/pow(c*mach*mach-1.0, d);
+ }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+// Based on the formulas in the US Air Force Aircraft Performance Flight Testing
+// Manual (AFFTC-TIH-99-01). In particular sections 4.6 to 4.8.
+
+double FGAtmosphere::MachFromImpactPressure(double qc, double p) const
+{
+ constexpr double a = 2.0/(SHRatio-1.0);
+ constexpr double b = (SHRatio-1.0)/SHRatio;
+ constexpr double c = 2.0/b;
+ constexpr double d = 0.5*a;
+ const double coeff = pow(0.5*(SHRatio+1.0), -0.25*c)*pow(0.5*(SHRatio+1.0)/SHRatio, -0.5*d);
+
+ double A = qc / p + 1;
+ double M = sqrt(a*(pow(A, b) - 1.0)); // Equation (4.12)
+
+ if (M > 1.0)
+ for (unsigned int i = 0; i<10; i++)
+ M = coeff*sqrt(A*pow(1 - 1.0 / (c*M*M), d)); // Equation (4.17)
+
+ return M;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGAtmosphere::VcalibratedFromMach(double mach, double altitude) const
+{
+ double p = GetPressure(altitude);
+ double qc = PitotTotalPressure(mach, p) - p;
+ return StdDaySLsoundspeed * MachFromImpactPressure(qc, StdDaySLpressure);
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGAtmosphere::MachFromVcalibrated(double vcas, double altitude) const
+{
+ double p = GetPressure(altitude);
+ double qc = PitotTotalPressure(vcas / StdDaySLsoundspeed, StdDaySLpressure) - StdDaySLpressure;
+ return MachFromImpactPressure(qc, p);
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 void FGAtmosphere::bind(void)
 {
  PropertyManager->Tie("atmosphere/T-R", this, &FGAtmosphere::GetTemperature);

src/models/FGAtmosphere.h

@@ -206,13 +206,48 @@ class JSBSIM_API FGAtmosphere : public FGModel {
 
  virtual double GetPressureAltitude() const {return PressureAltitude;}
 
+ /** Compute the total pressure in front of the Pitot tube. It uses the
+ * Rayleigh formula for supersonic speeds (See "Introduction to Aerodynamics
+ * of a Compressible Fluid - H.W. Liepmann, A.E. Puckett - Wiley & sons
+ * (1947)" §5.4 pp 75-80)
+ * @param mach The Mach number
+ * @param p Pressure in psf
+ * @return The total pressure in front of the Pitot tube in psf */
+ double PitotTotalPressure(double mach, double p) const;
+
+ /** Compute the Mach number from the differential pressure (qc) and the
+ * static pressure. Based on the formulas in the US Air Force Aircraft
+ * Performance Flight Testing Manual (AFFTC-TIH-99-01).
+ * @param qc The differential/impact pressure
+ * @param p Pressure in psf
+ * @return The Mach number */
+ double MachFromImpactPressure(double qc, double p) const;
+
+ /** Calculate the calibrated airspeed from the Mach number. Based on the
+ * formulas in the US Air Force Aircraft Performance Flight Testing
+ * Manual (AFFTC-TIH-99-01).
+ * @param mach The Mach number
+ * @param altitude The altitude above sea level (ASL) in feet.
+ * @return The calibrated airspeed (CAS) in ft/s
+ * */
+ double VcalibratedFromMach(double mach, double altitude) const;
+
+ /** Calculate the Mach number from the calibrated airspeed.Based on the
+ * formulas in the US Air Force Aircraft Performance Flight Testing
+ * Manual (AFFTC-TIH-99-01).
+ * @param vcas The calibrated airspeed (CAS) in ft/s
+ * @param altitude The altitude above sea level (ASL) in feet.
+ * @return The Mach number
+ * */
+ double MachFromVcalibrated(double vcas, double altitude) const;
+
  struct Inputs {
  double altitudeASL;
  } in;
 
  static constexpr double StdDaySLtemperature = 518.67;
  static constexpr double StdDaySLpressure = 2116.228;
- static const double StdDaySLsoundspeed;
+ const double StdDaySLsoundspeed;
 
 protected:
  // Sea level conditions

src/models/FGAuxiliary.cpp

@@ -47,6 +47,7 @@ INCLUDES
 #include "FGFDMExec.h"
 #include "input_output/FGPropertyManager.h"
 #include "FGInertial.h"
+#include "FGAtmosphere.h"
 
 using namespace std;
 
@@ -193,10 +194,10 @@ bool FGAuxiliary::Run(bool Holding)
  tat = in.Temperature*(1 + 0.2*Mach*Mach); // Total Temperature, isentropic flow
  tatc = RankineToCelsius(tat);
 
- pt = PitotTotalPressure(Mach, in.Pressure);
+ pt = FDMExec->GetAtmosphere()->PitotTotalPressure(Mach, in.Pressure);
 
  if (abs(Mach) > 0.0) {
- vcas = VcalibratedFromMach(Mach, in.Pressure);
+ vcas = FDMExec->GetAtmosphere()->VcalibratedFromMach(Mach, in.AltitudeASL);
  veas = sqrt(2 * qbar / in.DensitySL);
  }
  else

src/models/FGAuxiliary.h

@@ -260,6 +260,7 @@ class JSBSIM_API FGAuxiliary : public FGModel {
  double SoundSpeed;
  double KinematicViscosity;
  double DistanceAGL;
+ double AltitudeASL;
  double Wingspan;
  double Wingchord;
  double StandardGravity;