Move the Mach/VCAS conversion methods to FGAuxiliary (JSBSim-Team#898)

src/FGFDMExec.cpp

@@ -390,7 +390,6 @@ 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();
@@ -534,11 +533,12 @@ void FGFDMExec::LoadInputs(unsigned int idx)
 
 void FGFDMExec::LoadPlanetConstants(void)
 {
- Propagate->in.vOmegaPlanet = Inertial->GetOmegaPlanet();
- Accelerations->in.vOmegaPlanet = Inertial->GetOmegaPlanet();
- Propagate->in.SemiMajor = Inertial->GetSemimajor();
- Propagate->in.SemiMinor = Inertial->GetSemiminor();
- Auxiliary->in.StandardGravity = Inertial->GetStandardGravity();
+ Propagate->in.vOmegaPlanet = Inertial->GetOmegaPlanet();
+ Accelerations->in.vOmegaPlanet = Inertial->GetOmegaPlanet();
+ Propagate->in.SemiMajor = Inertial->GetSemimajor();
+ Propagate->in.SemiMinor = Inertial->GetSemiminor();
+ Auxiliary->in.StandardGravity = Inertial->GetStandardGravity();
+ Auxiliary->in.StdDaySLsoundspeed = Atmosphere->StdDaySLsoundspeed;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/initialization/FGInitialCondition.cpp

@@ -48,6 +48,7 @@ INCLUDES
 #include "models/FGInertial.h"
 #include "models/FGAtmosphere.h"
 #include "models/FGAccelerations.h"
+#include "models/FGAuxiliary.h"
 #include "input_output/FGXMLFileRead.h"
 #include "FGTrim.h"
 
@@ -64,6 +65,7 @@ FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec) : fdmex(FDMExec)
  if(FDMExec) {
  Atmosphere=fdmex->GetAtmosphere();
  Aircraft=fdmex->GetAircraft();
+ Auxiliary=fdmex->GetAuxiliary();
  } else {
  cout << "FGInitialCondition: This class requires a pointer to a valid FGFDMExec object" << endl;
  }
@@ -173,7 +175,8 @@ void FGInitialCondition::SetMachIC(double mach)
 void FGInitialCondition::SetVcalibratedKtsIC(double vcas)
 {
  double altitudeASL = GetAltitudeASLFtIC();
- double mach = Atmosphere->MachFromVcalibrated(fabs(vcas)*ktstofps, altitudeASL);
+ double pressure = Atmosphere->GetPressure(altitudeASL);
+ double mach = Auxiliary->MachFromVcalibrated(fabs(vcas)*ktstofps, pressure);
  double soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
 
  SetVtrueFpsIC(mach * soundSpeed);
@@ -675,12 +678,13 @@ 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 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL);
+ double vc0 = Auxiliary->VcalibratedFromMach(mach0, pressure);
  double ve0 = vt * sqrt(rho/rhoSL);
 
  switch(lastLatitudeSet) {
@@ -718,10 +722,11 @@ void FGInitialCondition::SetAltitudeAGLFtIC(double agl)
  altitudeASL = GetAltitudeASLFtIC();
  soundSpeed = Atmosphere->GetSoundSpeed(altitudeASL);
  rho = Atmosphere->GetDensity(altitudeASL);
+ pressure = Atmosphere->GetPressure(altitudeASL);
 
  switch(lastSpeedSet) {
  case setvc:
- mach0 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL);
+ mach0 = Auxiliary->MachFromVcalibrated(vc0, pressure);
  SetVtrueFpsIC(mach0 * soundSpeed);
  break;
  case setmach:
@@ -745,12 +750,13 @@ 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 = Atmosphere->VcalibratedFromMach(mach0, altitudeASL);
+ double vc0 = Auxiliary->VcalibratedFromMach(mach0, pressure);
  double ve0 = vt * sqrt(rho/rhoSL);
 
  switch(lastLatitudeSet) {
@@ -820,10 +826,11 @@ 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 = Atmosphere->MachFromVcalibrated(vc0, altitudeASL);
+ mach0 = Auxiliary->MachFromVcalibrated(vc0, pressure);
  SetVtrueFpsIC(mach0 * soundSpeed);
  break;
  case setmach:
@@ -953,10 +960,11 @@ 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 * Atmosphere->VcalibratedFromMach(mach, altitudeASL);
+ return fpstokts * Auxiliary->VcalibratedFromMach(mach, pressure);
 }
 
 //******************************************************************************

src/initialization/FGInitialCondition.h

@@ -62,6 +62,7 @@ class FGMatrix33;
 class FGColumnVector3;
 class FGAtmosphere;
 class FGAircraft;
+class FGAuxiliary;
 class FGPropertyManager;
 class Element;
 
@@ -705,6 +706,7 @@ class JSBSIM_API FGInitialCondition : public FGJSBBase
  FGFDMExec *fdmex;
  FGAtmosphere* Atmosphere;
  FGAircraft* Aircraft;
+ FGAuxiliary* Auxiliary;
 
  bool Load_v1(Element* document);
  bool Load_v2(Element* document);

src/models/FGAtmosphere.cpp

@@ -309,79 +309,6 @@ 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,48 +206,14 @@ 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;
  const double StdDaySLsoundspeed;
+ static constexpr double SHRatio = 1.4;
 
 protected:
  // Sea level conditions
@@ -323,8 +289,6 @@ class JSBSIM_API FGAtmosphere : public FGModel {
  static constexpr double Reng0 = Rstar / Mair;
  //@}
 
- static constexpr double SHRatio = 1.4;
-
  double Reng = Reng0;
 
  virtual void bind(void);

src/models/FGAuxiliary.cpp

@@ -61,9 +61,9 @@ CLASS IMPLEMENTATION
 FGAuxiliary::FGAuxiliary(FGFDMExec* fdmex) : FGModel(fdmex)
 {
  Name = "FGAuxiliary";
- pt = 2116.23; // ISA SL pressure
- tatc = 15.0; // ISA SL temperature
- tat = 518.67;
+ pt = FGAtmosphere::StdDaySLpressure;  // ISA SL pressure
+ tat = FGAtmosphere::StdDaySLtemperature; // ISA SL temperature
+ tatc = RankineToCelsius(tat);
 
  vcas = veas = 0.0;
  qbar = qbarUW = qbarUV = 0.0;
@@ -194,10 +194,10 @@ bool FGAuxiliary::Run(bool Holding)
  tat = in.Temperature*(1 + 0.2*Mach*Mach); // Total Temperature, isentropic flow
  tatc = RankineToCelsius(tat);
 
- pt = FDMExec->GetAtmosphere()->PitotTotalPressure(Mach, in.Pressure);
+ pt = PitotTotalPressure(Mach, in.Pressure);
 
  if (abs(Mach) > 0.0) {
- vcas = FDMExec->GetAtmosphere()->VcalibratedFromMach(Mach, in.AltitudeASL);
+ vcas = VcalibratedFromMach(Mach, in.Pressure);
  veas = sqrt(2 * qbar / in.DensitySL);
  }
  else
@@ -229,6 +229,82 @@ bool FGAuxiliary::Run(bool Holding)
  return false;
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGAuxiliary::PitotTotalPressure(double mach, double pressure) const
+{
+ constexpr double SHRatio = FGAtmosphere::SHRatio;
+ 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 pressure;
+ if (mach < 1) //calculate total pressure assuming isentropic flow
+ return pressure*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 pressure*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 FGAuxiliary::MachFromImpactPressure(double qc, double pressure) const
+{
+ constexpr double SHRatio = FGAtmosphere::SHRatio;
+ 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 / pressure + 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 FGAuxiliary::VcalibratedFromMach(double mach, double pressure) const
+{
+ double qc = PitotTotalPressure(mach, pressure) - pressure;
+ return in.StdDaySLsoundspeed * MachFromImpactPressure(qc, FGAtmosphere::StdDaySLpressure);
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGAuxiliary::MachFromVcalibrated(double vcas, double pressure) const
+{
+ constexpr double StdDaySLpressure = FGAtmosphere::StdDaySLpressure;
+ double qc = PitotTotalPressure(vcas / in.StdDaySLsoundspeed, StdDaySLpressure) - StdDaySLpressure;
+ return MachFromImpactPressure(qc, pressure);
+}
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 //
 // From Stevens and Lewis, "Aircraft Control and Simulation", 3rd Ed., the