namespace JSBSim {
-static const char *IdSrc = "$Id: FGJSBBase.cpp,v 1.29 2010/03/18 13:19:21 jberndt Exp $";
+static const char *IdSrc = "$Id: FGJSBBase.cpp,v 1.32 2011/10/22 14:38:30 bcoconni Exp $";
static const char *IdHdr = ID_JSBBASE;
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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';
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+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