Merge branch 'next' of git://gitorious.org/fg/flightgear into next

[flightgear.git] / src / FDM / JSBSim / models / FGAtmosphere.cpp

diff --git a/src/FDM/JSBSim/models/FGAtmosphere.cpp b/src/FDM/JSBSim/models/FGAtmosphere.cpp
index d2c6a3dac1adee4f3d40882e790c962da3c49388..990adb882818692fd8a11d2067d2c761f6d5b60b 100644 (file)
--- a/src/FDM/JSBSim/models/FGAtmosphere.cpp
+++ b/src/FDM/JSBSim/models/FGAtmosphere.cpp
@@ -61,7 +61,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGAtmosphere.cpp,v 1.40 2010/11/18 12:38:06 jberndt Exp $";
+static const char *IdSrc = "$Id: FGAtmosphere.cpp,v 1.45 2011/05/20 03:18:36 jberndt Exp $";
 static const char *IdHdr = ID_ATMOSPHERE;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -84,8 +84,7 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
   htab[7]=278385.0; //ft.
 
   MagnitudedAccelDt = MagnitudeAccel = Magnitude = 0.0;
-//  SetTurbType( ttCulp );
-  SetTurbType( ttNone );
+  SetTurbType( ttMilspec );
   TurbGain = 1.0;
   TurbRate = 10.0;
   Rhythmicity = 0.1;
@@ -124,6 +123,7 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
 
 FGAtmosphere::~FGAtmosphere()
 {
+  delete(POE_Table);
   Debug(1);
 }
 
@@ -131,8 +131,6 @@ FGAtmosphere::~FGAtmosphere()
 
 bool FGAtmosphere::InitModel(void)
 {
-  if (!FGModel::InitModel()) return false;
-
   UseInternal();  // this is the default
 
   Calculate(h);
@@ -150,10 +148,10 @@ bool FGAtmosphere::InitModel(void)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-bool FGAtmosphere::Run(void)
+bool FGAtmosphere::Run(bool Holding)
 {
-  if (FGModel::Run()) return true;
-  if (FDMExec->Holding()) return false;
+  if (FGModel::Run(Holding)) return true;
+  if (Holding) return false;
 
   RunPreFunctions();
 
@@ -482,45 +480,7 @@ void FGAtmosphere::Turbulence(void)
 
     break;
   }
-  case ttBerndt: { // This is very experimental and incomplete at the moment.
-
-    vDirectiondAccelDt(eX) = GaussianRandomNumber();
-    vDirectiondAccelDt(eY) = GaussianRandomNumber();
-    vDirectiondAccelDt(eZ) = GaussianRandomNumber();
-/*
-    MagnitudedAccelDt = GaussianRandomNumber();
-    MagnitudeAccel    += MagnitudedAccelDt * DeltaT;
-    Magnitude         += MagnitudeAccel * DeltaT;
-*/
-    Magnitude         += GaussianRandomNumber() * DeltaT;
-
-    vDirectiondAccelDt.Normalize();
-    vDirectionAccel += TurbRate * vDirectiondAccelDt * DeltaT;
-    vDirectionAccel.Normalize();
-    vDirection      += vDirectionAccel*DeltaT;
-
-    // Diminish z-vector within two wingspans of the ground
-    if (HOverBMAC < 2.0) vDirection(eZ) *= HOverBMAC / 2.0;
-
-    vDirection.Normalize();
-
-    vTurbulenceNED = TurbGain*Magnitude * vDirection;
-    vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
 
-    vBodyTurbGrad = Tl2b * vTurbulenceGrad;
-    vTurbPQR(eP) = vBodyTurbGrad(eY) / wingspan;
-    if (HTailArm > 0)
-      vTurbPQR(eQ) = vBodyTurbGrad(eZ) / HTailArm;
-    else
-      vTurbPQR(eQ) = vBodyTurbGrad(eZ) / 10.0;
-
-    if (VTailArm > 0)
-      vTurbPQR(eR) = vBodyTurbGrad(eX) / VTailArm;
-    else
-      vTurbPQR(eR) = vBodyTurbGrad(eX)/10.0;
-
-    break;
-  }
   case ttCulp: { 
 
     vTurbPQR(eP) = wind_from_clockwise;
@@ -572,8 +532,11 @@ void FGAtmosphere::Turbulence(void)
   }
   case ttMilspec:
   case ttTustin: {
+    double V = FDMExec->GetAuxiliary()->GetVt(); // true airspeed in ft/s
+
     // an index of zero means turbulence is disabled
-    if (probability_of_exceedence_index == 0) {
+    // airspeed occurs as divisor in the code below
+    if (probability_of_exceedence_index == 0 || V == 0) {
       vTurbulenceNED(1) = vTurbulenceNED(2) = vTurbulenceNED(3) = 0.0;
       vTurbPQR(1) = vTurbPQR(2) = vTurbPQR(3) = 0.0;
       return;
@@ -582,10 +545,11 @@ void FGAtmosphere::Turbulence(void)
     // Turbulence model according to MIL-F-8785C (Flying Qualities of Piloted Aircraft)
     double
       h = FDMExec->GetPropagate()->GetDistanceAGL(),
-      V = FDMExec->GetAuxiliary()->GetVt(), // true airspeed in ft/s
       b_w = wingspan,
       L_u, L_w, sig_u, sig_w;
 
+      if (b_w == 0.) b_w = 30.;
+
     // clip height functions at 10 ft
     if (h <= 10.) h = 10;
 
@@ -644,6 +608,10 @@ void FGAtmosphere::Turbulence(void)
         C_BLq = 1/tau_q/tan(T_V/2/tau_q), // eq. (24)
         C_BLr = 1/tau_r/tan(T_V/2/tau_r); // eq. (26)
 
+      // all values calculated so far are strictly positive, except for
+      // the random numbers nu_*. This means that in the code below, all
+      // divisors are strictly positive, too, and no floating point
+      // exception should occur.
       xi_u = -(1 - C_BL*tau_u)/(1 + C_BL*tau_u)*xi_u_km1
            + sig_u*sqrt(2*tau_u/T_V)/(1 + C_BL*tau_u)*(nu_u + nu_u_km1); // eq. (18)
       xi_v = -2*(sqr(omega_v) - sqr(C_BL))/sqr(omega_v + C_BL)*xi_v_km1