- switch (turbType) {
- case ttStandard: {
- vDirectiondAccelDt(eX) = 1 - 2.0*(double(rand())/double(RAND_MAX));
- vDirectiondAccelDt(eY) = 1 - 2.0*(double(rand())/double(RAND_MAX));
- vDirectiondAccelDt(eZ) = 1 - 2.0*(double(rand())/double(RAND_MAX));
-
- MagnitudedAccelDt = 1 - 2.0*(double(rand())/double(RAND_MAX)) - Magnitude;
- // Scale the magnitude so that it moves
- // away from the peaks
- MagnitudedAccelDt = ((MagnitudedAccelDt - Magnitude) /
- (1 + fabs(Magnitude)));
- MagnitudeAccel += MagnitudedAccelDt*rate*TurbRate*State->Getdt();
- Magnitude += MagnitudeAccel*rate*State->Getdt();
-
- vDirectiondAccelDt.Normalize();
-
- // deemphasise non-vertical forces
- vDirectiondAccelDt(eX) = square_signed(vDirectiondAccelDt(eX));
- vDirectiondAccelDt(eY) = square_signed(vDirectiondAccelDt(eY));
-
- vDirectionAccel += vDirectiondAccelDt*rate*TurbRate*State->Getdt();
- vDirectionAccel.Normalize();
- vDirection += vDirectionAccel*rate*State->Getdt();
-
- vDirection.Normalize();
-
- // Diminish turbulence within three wingspans
- // of the ground
- vTurbulence = TurbGain * Magnitude * vDirection;
- double HOverBMAC = Auxiliary->GetHOverBMAC();
- if (HOverBMAC < 3.0)
- vTurbulence *= (HOverBMAC / 3.0) * (HOverBMAC / 3.0);
-
- vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
-
- vBodyTurbGrad = Propagate->GetTl2b()*vTurbulenceGrad;
-
- if (Aircraft->GetWingSpan() > 0) {
- vTurbPQR(eP) = vBodyTurbGrad(eY)/Aircraft->GetWingSpan();
- } else {
- vTurbPQR(eP) = vBodyTurbGrad(eY)/30.0;
- }
-// if (Aircraft->GetHTailArm() != 0.0)
-// vTurbPQR(eQ) = vBodyTurbGrad(eZ)/Aircraft->GetHTailArm();
-// else
-// vTurbPQR(eQ) = vBodyTurbGrad(eZ)/10.0;
-
- if (Aircraft->GetVTailArm() > 0)
- vTurbPQR(eR) = vBodyTurbGrad(eX)/Aircraft->GetVTailArm();
- else
- vTurbPQR(eR) = vBodyTurbGrad(eX)/10.0;
-
- // Clear the horizontal forces
- // actually felt by the plane, now
- // that we've used them to calculate
- // moments.
- vTurbulence(eX) = 0.0;
- vTurbulence(eY) = 0.0;