Date started: 07/01/01
Purpose: Encapsulates the JSBBase object
- ------------- Copyright (C) 2001 Jon S. Berndt (jsb@hal-pc.org) -------------
+ ------------- Copyright (C) 2001 Jon S. Berndt (jon@jsbsim.org) -------------
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the Free Software
#include "FGJSBBase.h"
#include <iostream>
+#include <sstream>
+#include <cstdlib>
namespace JSBSim {
-static const char *IdSrc = "$Id$";
+static const char *IdSrc = "$Id: FGJSBBase.cpp,v 1.32 2011/10/22 14:38:30 bcoconni Exp $";
static const char *IdHdr = ID_JSBBASE;
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
char FGJSBBase::fgdef[6] = {'\0' };
#endif
-const double FGJSBBase::radtodeg = 57.29578;
-const double FGJSBBase::degtorad = 1.745329E-2;
+const double FGJSBBase::radtodeg = 57.295779513082320876798154814105;
+const double FGJSBBase::degtorad = 0.017453292519943295769236907684886;
const double FGJSBBase::hptoftlbssec = 550.0;
const double FGJSBBase::psftoinhg = 0.014138;
const double FGJSBBase::psftopa = 47.88;
const double FGJSBBase::m3toft3 = 1.0/(fttom*fttom*fttom);
const double FGJSBBase::inhgtopa = 3386.38;
const double FGJSBBase::fttom = 0.3048;
-double FGJSBBase::Reng = 1716.0;
+double FGJSBBase::Reng = 1716.56; // Gas constant for Air (ft-lb/slug-R)
+double FGJSBBase::Rstar = 1545.348; // Universal gas constant
+double FGJSBBase::Mair = 28.9645; //
const double FGJSBBase::SHRatio = 1.40;
// Note that definition of lbtoslug by the inverse of slugtolb and not
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-void FGJSBBase::disableHighLighting(void) {
+void FGJSBBase::disableHighLighting(void)
+{
highint[0]='\0';
halfint[0]='\0';
normint[0]='\0';
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+string FGJSBBase::CreateIndexedPropertyName(const string& Property, int index)
+{
+ std::ostringstream buf;
+ buf << Property << '[' << index << ']';
+ return buf.str();
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGJSBBase::GaussianRandomNumber(void)
+{
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+
+ if (phase == 0) {
+ V1 = V2 = S = X = 0.0;
+
+ do {
+ double U1 = (double)rand() / RAND_MAX;
+ double U2 = (double)rand() / RAND_MAX;
+
+ V1 = 2 * U1 - 1;
+ V2 = 2 * U2 - 1;
+ S = V1 * V1 + V2 * V2;
+ } while(S >= 1 || S == 0);
+
+ X = V1 * sqrt(-2 * log(S) / S);
+ } else
+ X = V2 * sqrt(-2 * log(S) / S);
+
+ phase = 1 - phase;
+
+ return X;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGJSBBase::VcalibratedFromMach(double mach, double p, double psl, double rhosl)
+{
+ double pt,A;
+
+ if (mach < 0) mach=0;
+ if (mach < 1) //calculate total pressure assuming isentropic flow
+ pt=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
+
+ pt = p*166.92158*pow(mach,7.0)/pow(7*mach*mach-1,2.5);
+ }
+
+ A = pow(((pt-p)/psl+1),0.28571);
+ return sqrt(7*psl/rhosl*(A-1));
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGJSBBase::MachFromVcalibrated(double vcas, double p, double psl, double rhosl)
+{
+ double pt = p + psl*(pow(1+vcas*vcas*rhosl/(7.0*psl),3.5)-1);
+
+ if (pt/p < 1.89293)
+ return sqrt(5.0*(pow(pt/p, 0.2857143) -1)); // Mach < 1
+ else {
+ // Mach >= 1
+ double mach = sqrt(0.77666*pt/p); // Initial guess is based on a quadratic approximation of the Rayleigh formula
+ double delta = 1.;
+ double target = pt/(166.92158*p);
+ int iter = 0;
+
+ // Find the root with Newton-Raphson. Since the differential is never zero,
+ // the function is monotonic and has only one root with a multiplicity of one.
+ // Convergence is certain.
+ while (delta > 1E-5 && iter < 10) {
+ double m2 = mach*mach; // Mach^2
+ double m6 = m2*m2*m2; // Mach^6
+ delta = mach*m6/pow(7.0*m2-1.0,2.5) - target;
+ double diff = 7.0*m6*(2.0*m2-1)/pow(7.0*m2-1.0,3.5); // Never zero when Mach >= 1
+ mach -= delta/diff;
+ iter++;
+ }
+
+ return mach;
+ }
+}
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
} // namespace JSBSim