45#include "FGAuxiliary.h" 
   46#include "initialization/FGInitialCondition.h" 
   48#include "input_output/FGPropertyManager.h" 
   49#include "FGInertial.h" 
   50#include "FGAtmosphere.h" 
   64  pt = FGAtmosphere::StdDaySLpressure;     
 
   65  tat = FGAtmosphere::StdDaySLtemperature; 
 
   69  qbar = qbarUW = qbarUV = 0.0;
 
   73  gamma = Vt = Vground = 0.0;
 
   75  hoverbmac = hoverbcg = 0.0;
 
   79  vPilotAccel.InitMatrix();
 
   80  vPilotAccelN.InitMatrix();
 
   81  vAeroUVW.InitMatrix();
 
   82  vAeroPQR.InitMatrix();
 
   83  vMachUVW.InitMatrix();
 
   84  vEulerRates.InitMatrix();
 
 
   93bool FGAuxiliary::InitModel(
void)
 
   95  if (!FGModel::InitModel()) 
return false;
 
  102  qbar = qbarUW = qbarUV = 0.0;
 
  106  gamma = Vt = Vground = 0.0;
 
  108  hoverbmac = hoverbcg = 0.0;
 
  112  vPilotAccel.InitMatrix();
 
  113  vPilotAccelN.InitMatrix();
 
  114  vAeroUVW.InitMatrix();
 
  115  vAeroPQR.InitMatrix();
 
  116  vMachUVW.InitMatrix();
 
  117  vEulerRates.InitMatrix();
 
  134  if (Holding) 
return false;
 
  138  vEulerRates(eTht) = in.vPQR(eQ)*in.CosPhi - in.vPQR(eR)*in.SinPhi;
 
  139  if (in.CosTht != 0.0) {
 
  140    vEulerRates(ePsi) = (in.vPQR(eQ)*in.SinPhi + in.vPQR(eR)*in.CosPhi)/in.CosTht;
 
  141    vEulerRates(ePhi) = in.vPQR(eP) + vEulerRates(ePsi)*in.SinTht;
 
  145  vAeroPQR = in.vPQR - in.TurbPQR;
 
  146  vAeroUVW = in.vUVW - in.Tl2b * in.TotalWindNED;
 
  148  alpha = beta = adot = bdot = 0;
 
  149  double AeroU2 = vAeroUVW(eU)*vAeroUVW(eU);
 
  150  double AeroV2 = vAeroUVW(eV)*vAeroUVW(eV);
 
  151  double AeroW2 = vAeroUVW(eW)*vAeroUVW(eW);
 
  152  double mUW = AeroU2 + AeroW2;
 
  154  double Vt2 = mUW + AeroV2;
 
  158    beta = atan2(vAeroUVW(eV), sqrt(mUW));
 
  161      alpha = atan2(vAeroUVW(eW), vAeroUVW(eU));
 
  162      double Vtdot = (vAeroUVW(eU)*in.vUVWdot(eU) + vAeroUVW(eV)*in.vUVWdot(eV) + vAeroUVW(eW)*in.vUVWdot(eW))/Vt;
 
  163      adot = (vAeroUVW(eU)*in.vUVWdot(eW) - vAeroUVW(eW)*in.vUVWdot(eU))/mUW;
 
  164      bdot = (in.vUVWdot(eV)*Vt - vAeroUVW(eV)*Vtdot)/(Vt*sqrt(mUW));
 
  168  UpdateWindMatrices();
 
  170  Re = Vt * in.Wingchord / in.KinematicViscosity;
 
  172  double densityD2 = 0.5*in.Density;
 
  174  qbar = densityD2 * Vt2;
 
  175  qbarUW = densityD2 * (mUW);
 
  176  qbarUV = densityD2 * (AeroU2 + AeroV2);
 
  177  Mach = Vt / in.SoundSpeed;
 
  178  MachU = vMachUVW(eU) = vAeroUVW(eU) / in.SoundSpeed;
 
  179  vMachUVW(eV) = vAeroUVW(eV) / in.SoundSpeed;
 
  180  vMachUVW(eW) = vAeroUVW(eW) / in.SoundSpeed;
 
  184  Vground = sqrt( in.vVel(eNorth)*in.vVel(eNorth) + in.vVel(eEast)*in.vVel(eEast) );
 
  186  psigt = atan2(in.vVel(eEast), in.vVel(eNorth));
 
  187  if (psigt < 0.0) psigt += 2*M_PI;
 
  188  gamma = atan2(-in.vVel(eDown), Vground);
 
  190  tat = in.Temperature*(1 + 0.2*Mach*Mach); 
 
  195  if (abs(Mach) > 0.0) {
 
  197    veas = sqrt(2 * qbar / FGAtmosphere::StdDaySLdensity);
 
  202  vPilotAccel.InitMatrix();
 
  203  vNcg = in.vBodyAccel/in.StandardGravity;
 
  208  vPilotAccel = in.vBodyAccel + in.vPQRidot * in.ToEyePt;
 
  209  vPilotAccel += in.vPQRi * (in.vPQRi * in.ToEyePt);
 
  211  vNwcg = mTb2w * vNcg;
 
  212  vNwcg(eZ) = 1.0 - vNwcg(eZ);
 
  214  vPilotAccelN = vPilotAccel / in.StandardGravity;
 
  220  hoverbcg = in.DistanceAGL / in.Wingspan;
 
  223  hoverbmac = (in.DistanceAGL - vMac(3)) / in.Wingspan;
 
 
  232  constexpr double SHRatio = FGAtmosphere::SHRatio;
 
  233  constexpr double a = (SHRatio-1.0) / 2.0;
 
  234  constexpr double b = SHRatio / (SHRatio-1.0);
 
  235  constexpr double c = 2.0*b;
 
  236  constexpr double d = 1.0 / (SHRatio-1.0);
 
  237  const double coeff = pow(0.5*(SHRatio+1.0), b)
 
  238                     * pow((SHRatio+1.0)/(SHRatio-1.0), d);
 
  240  if (mach < 0) 
return pressure;
 
  242    return pressure*pow((1.0 + a*mach*mach), b);
 
  258    return pressure*coeff*pow(mach, c)/pow(c*mach*mach-1.0, d);
 
 
  269  constexpr double SHRatio = FGAtmosphere::SHRatio;
 
  270  constexpr double a = 2.0/(SHRatio-1.0);
 
  271  constexpr double b = (SHRatio-1.0)/SHRatio;
 
  272  constexpr double c = 2.0/b;
 
  273  constexpr double d = 0.5*a;
 
  274  const double coeff = pow(0.5*(SHRatio+1.0), -0.25*c)
 
  275                     * pow(0.5*(SHRatio+1.0)/SHRatio, -0.5*d);
 
  277  double A = qc / pressure + 1;
 
  278  double M = sqrt(a*(pow(A, b) - 1.0));  
 
  281    for (
unsigned int i = 0; i<10; i++)
 
  282      M = coeff*sqrt(A*pow(1 - 1.0 / (c*M*M), d));  
 
 
  299  constexpr double StdDaySLpressure = FGAtmosphere::StdDaySLpressure;
 
  300  double qc = 
PitotTotalPressure(vcas / in.StdDaySLsoundspeed, StdDaySLpressure) - StdDaySLpressure;
 
 
  320void FGAuxiliary::UpdateWindMatrices(
void)
 
  322  double ca, cb, sa, sb;
 
  344double FGAuxiliary::GetNlf(
void)
 const 
  347    return (in.vFw(3))/(in.Mass*
slugtolb);
 
  354double FGAuxiliary::GetLongitudeRelativePosition(
void)
 const 
  362double FGAuxiliary::GetLatitudeRelativePosition(
void)
 const 
  365                                    FDMExec->
GetIC()->GetGeodLatitudeRadIC())* fttom;
 
  370double FGAuxiliary::GetDistanceRelativePosition(
void)
 const 
  372  auto ic = FDMExec->
GetIC();
 
  374                                    ic->GetGeodLatitudeRadIC())* fttom;
 
  379void FGAuxiliary::bind(
void)
 
  381  typedef double (FGAuxiliary::*PMF)(int) const;
 
  382  typedef double (FGAuxiliary::*PF)(void) const;
 
  384  PropertyManager->Tie(
"propulsion/tat-c", 
this, &FGAuxiliary::GetTAT_C);
 
  393  PropertyManager->Tie(
"velocities/p-aero-rad_sec", 
this, eX, (PMF)&FGAuxiliary::GetAeroPQR);
 
  394  PropertyManager->Tie(
"velocities/q-aero-rad_sec", 
this, eY, (PMF)&FGAuxiliary::GetAeroPQR);
 
  395  PropertyManager->Tie(
"velocities/r-aero-rad_sec", 
this, eZ, (PMF)&FGAuxiliary::GetAeroPQR);
 
  396  PropertyManager->Tie(
"velocities/phidot-rad_sec", 
this, ePhi, (PMF)&FGAuxiliary::GetEulerRates);
 
  397  PropertyManager->Tie(
"velocities/thetadot-rad_sec", 
this, eTht, (PMF)&FGAuxiliary::GetEulerRates);
 
  398  PropertyManager->Tie(
"velocities/psidot-rad_sec", 
this, ePsi, (PMF)&FGAuxiliary::GetEulerRates);
 
  399  PropertyManager->Tie(
"velocities/u-aero-fps", 
this, eU, (PMF)&FGAuxiliary::GetAeroUVW);
 
  400  PropertyManager->Tie(
"velocities/v-aero-fps", 
this, eV, (PMF)&FGAuxiliary::GetAeroUVW);
 
  401  PropertyManager->Tie(
"velocities/w-aero-fps", 
this, eW, (PMF)&FGAuxiliary::GetAeroUVW);
 
  405  PropertyManager->Tie(
"accelerations/a-pilot-x-ft_sec2", 
this, eX, (PMF)&FGAuxiliary::GetPilotAccel);
 
  406  PropertyManager->Tie(
"accelerations/a-pilot-y-ft_sec2", 
this, eY, (PMF)&FGAuxiliary::GetPilotAccel);
 
  407  PropertyManager->Tie(
"accelerations/a-pilot-z-ft_sec2", 
this, eZ, (PMF)&FGAuxiliary::GetPilotAccel);
 
  408  PropertyManager->Tie(
"accelerations/n-pilot-x-norm", 
this, eX, (PMF)&FGAuxiliary::GetNpilot);
 
  409  PropertyManager->Tie(
"accelerations/n-pilot-y-norm", 
this, eY, (PMF)&FGAuxiliary::GetNpilot);
 
  410  PropertyManager->Tie(
"accelerations/n-pilot-z-norm", 
this, eZ, (PMF)&FGAuxiliary::GetNpilot);
 
  414  PropertyManager->Tie(
"forces/load-factor", 
this, &FGAuxiliary::GetNlf);
 
  415  PropertyManager->Tie(
"aero/alpha-rad", 
this, (PF)&FGAuxiliary::Getalpha);
 
  416  PropertyManager->Tie(
"aero/beta-rad", 
this, (PF)&FGAuxiliary::Getbeta);
 
  417  PropertyManager->Tie(
"aero/mag-beta-rad", 
this, (PF)&FGAuxiliary::GetMagBeta);
 
  418  PropertyManager->Tie(
"aero/alpha-deg", 
this, inDegrees, (PMF)&FGAuxiliary::Getalpha);
 
  419  PropertyManager->Tie(
"aero/beta-deg", 
this, inDegrees, (PMF)&FGAuxiliary::Getbeta);
 
  420  PropertyManager->Tie(
"aero/mag-beta-deg", 
this, inDegrees, (PMF)&FGAuxiliary::GetMagBeta);
 
  421  PropertyManager->Tie(
"aero/Re", 
this, &FGAuxiliary::GetReynoldsNumber);
 
  422  PropertyManager->Tie(
"aero/qbar-psf", 
this, &FGAuxiliary::Getqbar);
 
  423  PropertyManager->Tie(
"aero/qbarUW-psf", 
this, &FGAuxiliary::GetqbarUW);
 
  424  PropertyManager->Tie(
"aero/qbarUV-psf", 
this, &FGAuxiliary::GetqbarUV);
 
  425  PropertyManager->Tie(
"aero/alphadot-rad_sec", 
this, (PF)&FGAuxiliary::Getadot);
 
  426  PropertyManager->Tie(
"aero/betadot-rad_sec", 
this, (PF)&FGAuxiliary::Getbdot);
 
  427  PropertyManager->Tie(
"aero/alphadot-deg_sec", 
this, inDegrees, (PMF)&FGAuxiliary::Getadot);
 
  428  PropertyManager->Tie(
"aero/betadot-deg_sec", 
this, inDegrees, (PMF)&FGAuxiliary::Getbdot);
 
  429  PropertyManager->Tie(
"aero/h_b-cg-ft", 
this, &FGAuxiliary::GetHOverBCG);
 
  430  PropertyManager->Tie(
"aero/h_b-mac-ft", 
this, &FGAuxiliary::GetHOverBMAC);
 
  431  PropertyManager->Tie(
"flight-path/gamma-rad", 
this, &FGAuxiliary::GetGamma);
 
  432  PropertyManager->Tie(
"flight-path/gamma-deg", 
this, inDegrees, (PMF)&FGAuxiliary::GetGamma);
 
  433  PropertyManager->Tie(
"flight-path/psi-gt-rad", 
this, &FGAuxiliary::GetGroundTrack);
 
  435  PropertyManager->Tie(
"position/distance-from-start-lon-mt", 
this, &FGAuxiliary::GetLongitudeRelativePosition);
 
  436  PropertyManager->Tie(
"position/distance-from-start-lat-mt", 
this, &FGAuxiliary::GetLatitudeRelativePosition);
 
  437  PropertyManager->Tie(
"position/distance-from-start-mag-mt", 
this, &FGAuxiliary::GetDistanceRelativePosition);
 
  445double FGAuxiliary::BadUnits(
void)
 const 
  447  cerr << 
"Bad units" << endl; 
return 0.0;
 
  469void FGAuxiliary::Debug(
int from)
 
  471  if (debug_lvl <= 0) 
return;
 
  478  if (debug_lvl & 2 ) { 
 
  479    if (from == 0) cout << 
"Instantiated: FGAuxiliary" << endl;
 
  480    if (from == 1) cout << 
"Destroyed:    FGAuxiliary" << endl;
 
  482  if (debug_lvl & 4 ) { 
 
  484  if (debug_lvl & 8 ) { 
 
  486  if (debug_lvl & 16) { 
 
  487    if (Mach > 100 || Mach < 0.00)
 
  488      cout << 
"FGPropagate::Mach is out of bounds: " << Mach << endl;
 
  489    if (qbar > 1e6 || qbar < 0.00)
 
  490      cout << 
"FGPropagate::qbar is out of bounds: " << qbar << endl;
 
  492  if (debug_lvl & 64) {
 
double PitotTotalPressure(double mach, double pressure) const
Compute the total pressure in front of the Pitot tube.
 
double VcalibratedFromMach(double mach, double pressure) const
Calculate the calibrated airspeed from the Mach number.
 
double GetVtrueFPS() const
Returns the true airspeed in feet per second.
 
double GetVcalibratedFPS(void) const
Returns Calibrated airspeed in feet/second.
 
double GetVground(void) const
Gets the ground speed in feet per second.
 
double GetVcalibratedKTS(void) const
Returns Calibrated airspeed in knots.
 
double GetMach(void) const
Gets the Mach number.
 
double GetTotalPressure(void) const
Returns the total pressure.
 
double GetVtrueKTS() const
Returns the true airspeed in knots.
 
double GetNz(void) const
The vertical acceleration in g's of the aircraft center of gravity.
 
double GetVequivalentFPS(void) const
Returns equivalent airspeed in feet/second.
 
double GetMachU(void) const
The mach number calculated using the vehicle X axis velocity.
 
bool Run(bool Holding) override
Runs the Auxiliary routines; called by the Executive Can pass in a value indicating if the executive ...
 
double GetTotalTemperature(void) const
Returns the total temperature.
 
double GetVt(void) const
Gets the magnitude of total vehicle velocity including wind effects in feet per second.
 
double GetNy(void) const
The lateral acceleration in g's of the aircraft center of gravity.
 
double GetNx(void) const
The longitudinal acceleration in g's of the aircraft center of gravity.
 
~FGAuxiliary()
Destructor.
 
double MachFromVcalibrated(double vcas, double pressure) const
Calculate the Mach number from the calibrated airspeed.Based on the formulas in the US Air Force Airc...
 
double MachFromImpactPressure(double qc, double p) const
Compute the Mach number from the differential pressure (qc) and the static pressure.
 
FGAuxiliary(FGFDMExec *Executive)
Constructor.
 
double GetVequivalentKTS(void) const
Returns equivalent airspeed in knots.
 
This class implements a 3 element column vector.
 
Encapsulates the JSBSim simulation executive.
 
std::shared_ptr< FGInitialCondition > GetIC(void) const
Returns a pointer to the FGInitialCondition object.
 
static constexpr double slugtolb
Note that definition of lbtoslug by the inverse of slugtolb and not to a different constant you can a...
 
static constexpr double RankineToCelsius(double rankine)
Converts from degrees Rankine to degrees Celsius.
 
double GetGeodLatitudeRad(void) const
Get the GEODETIC latitude in radians.
 
double GetLongitude() const
Get the longitude.
 
double GetRadius() const
Get the distance from the center of the earth in feet.
 
FGLocation LocalToLocation(const FGColumnVector3 &lvec) const
Conversion from Local frame coordinates to a location in the earth centered and fixed frame.
 
double GetLatitudeDeg() const
Get the GEOCENTRIC latitude in degrees.
 
double GetLongitudeDeg() const
Get the longitude.
 
double GetDistanceTo(double target_longitude, double target_latitude) const
Get the geodetic distance between the current location and a given location.
 
FGMatrix33 Transposed(void) const
Transposed matrix.
 
Base class for all scheduled JSBSim models.
 
virtual bool Run(bool Holding)
Runs the model; called by the Executive.