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;
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
FGAtmosphere::~FGAtmosphere()
{
+ delete(POE_Table);
Debug(1);
}
bool FGAtmosphere::InitModel(void)
{
- if (!FGModel::InitModel()) return false;
-
UseInternal(); // this is the default
Calculate(h);
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-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();
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;
}
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;
// 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;
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