]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/JSBSim/models/FGAtmosphere.cpp
Sync. with JSBSim CVS
[flightgear.git] / src / FDM / JSBSim / models / FGAtmosphere.cpp
index 2441f843f7a1ac21a473b07c445e8fb7659b8723..88ffbd5e0fb8e8d06cb3b65c897ab77325ac7781 100644 (file)
  ------------- Copyright (C) 1999  Jon S. Berndt (jsb@hal-pc.org) -------------
 
  This program is free software; you can redistribute it and/or modify it under
- the terms of the GNU General Public License as published by the Free Software
+ the terms of the GNU Lesser General Public License as published by the Free Software
  Foundation; either version 2 of the License, or (at your option) any later
  version.
 
  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  details.
 
- You should have received a copy of the GNU General Public License along with
+ You should have received a copy of the GNU Lesser General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  Place - Suite 330, Boston, MA  02111-1307, USA.
 
- Further information about the GNU General Public License can also be found on
+ Further information about the GNU Lesser General Public License can also be found on
  the world wide web at http://www.gnu.org.
 
 FUNCTIONAL DESCRIPTION
@@ -64,7 +64,6 @@ static const char *IdHdr = ID_ATMOSPHERE;
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-
 FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
 {
   Name = "FGAtmosphere";
@@ -81,14 +80,18 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
   htab[7]=259186.352; //ft.
 
   MagnitudedAccelDt = MagnitudeAccel = Magnitude = 0.0;
-//   turbType = ttNone;
-  turbType = ttStandard;
-//   turbType = ttBerndt;
+  SetTurbType( ttCulp );
   TurbGain = 0.0;
-  TurbRate = 1.0;
+  TurbRate = 1.7;
+  Rhythmicity = 0.1;
+  spike = target_time = strength = 0.0;
+  wind_from_clockwise = 0.0;
 
   T_dev_sl = T_dev = delta_T = 0.0;
   StandardTempOnly = false;
+  first_pass = true;
+  vGustNED.InitMatrix();
+  vTurbulenceNED.InitMatrix();
 
   bind();
   Debug(0);
@@ -98,7 +101,6 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
 
 FGAtmosphere::~FGAtmosphere()
 {
-  unbind();
   Debug(1);
 }
 
@@ -106,7 +108,7 @@ FGAtmosphere::~FGAtmosphere()
 
 bool FGAtmosphere::InitModel(void)
 {
-  FGModel::InitModel();
+  if (!FGModel::InitModel()) return false;
 
   UseInternal();  // this is the default
 
@@ -151,9 +153,8 @@ bool FGAtmosphere::Run(void)
 void FGAtmosphere::Calculate(double altitude)
 {
   double slope, reftemp, refpress;
-  int i = 0;
+  int i = lastIndex;
 
-  i = lastIndex;
   if (altitude < htab[lastIndex]) {
     if (altitude <= 0) {
       i = 0;
@@ -255,17 +256,18 @@ void FGAtmosphere::Calculate(double altitude)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Calculate parameters derived from T, P and rho
+// Sum gust and turbulence values in NED frame into the wind vector.
 
 void FGAtmosphere::CalculateDerived(void)
 {
   T_dev = (*temperature) - GetTemperature(h);
   density_altitude = h + T_dev * 66.7;
 
-  if (turbType == ttStandard) {
-    Turbulence();
-    vWindNED += vTurbulence;
-  }
-  if (vWindNED(1) != 0.0) psiw = atan2( vWindNED(2), vWindNED(1) );
+  if (turbType == ttStandard || ttCulp) Turbulence();
+
+  vTotalWindNED = vWindNED + vGustNED + vTurbulenceNED;
+
+  if (vWindNED(eX) != 0.0) psiw = atan2( vWindNED(eY), vWindNED(eX) );
   if (psiw < 0) psiw += 2*M_PI;
 
   soundspeed = sqrt(SHRatio*Reng*(*temperature));
@@ -325,10 +327,42 @@ static inline double square_signed (double value)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+void FGAtmosphere::SetWindspeed(double speed)
+{
+  if (vWindNED.Magnitude() == 0.0) {
+    psiw = 0.0;
+    vWindNED(eNorth) = speed;
+  } else {
+    vWindNED(eNorth) = speed * cos(psiw);
+    vWindNED(eEast) = speed * sin(psiw);
+    vWindNED(eDown) = 0.0;
+  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGAtmosphere::GetWindspeed(void) const
+{
+  return vWindNED.Magnitude();
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void FGAtmosphere::SetWindPsi(double dir)
+{
+  double mag = GetWindspeed();
+  psiw = dir;
+  SetWindspeed(mag);  
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 void FGAtmosphere::Turbulence(void)
 {
   switch (turbType) {
   case ttStandard: {
+    TurbGain = TurbGain * TurbGain * 100.0;
+
     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));
@@ -340,6 +374,7 @@ void FGAtmosphere::Turbulence(void)
                          (1 + fabs(Magnitude)));
     MagnitudeAccel    += MagnitudedAccelDt*rate*TurbRate*State->Getdt();
     Magnitude         += MagnitudeAccel*rate*State->Getdt();
+    Magnitude          = fabs(Magnitude);
 
     vDirectiondAccelDt.Normalize();
 
@@ -355,13 +390,20 @@ void FGAtmosphere::Turbulence(void)
 
                                 // Diminish turbulence within three wingspans
                                 // of the ground
-    vTurbulence = TurbGain * Magnitude * vDirection;
+    vTurbulenceNED = TurbGain * Magnitude * vDirection;
     double HOverBMAC = Auxiliary->GetHOverBMAC();
     if (HOverBMAC < 3.0)
-        vTurbulence *= (HOverBMAC / 3.0) * (HOverBMAC / 3.0);
+        vTurbulenceNED *= (HOverBMAC / 3.0) * (HOverBMAC / 3.0);
 
-    vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
+    // I don't believe these next two statements calculate the proper gradient over
+    // the aircraft body. One reason is because this has no relationship with the
+    // orientation or velocity of the aircraft, which it must have. What is vTurbulenceGrad
+    // supposed to represent? And the direction and magnitude of the turbulence can change,
+    // so both accelerations need to be accounted for, no?
+
+    // Need to determine the turbulence change in body axes between two time points.
 
+    vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
     vBodyTurbGrad = Propagate->GetTl2b()*vTurbulenceGrad;
 
     if (Aircraft->GetWingSpan() > 0) {
@@ -383,12 +425,16 @@ void FGAtmosphere::Turbulence(void)
                                 // actually felt by the plane, now
                                 // that we've used them to calculate
                                 // moments.
-    vTurbulence(eX) = 0.0;
-    vTurbulence(eY) = 0.0;
+                                // Why? (JSB)
+//    vTurbulenceNED(eX) = 0.0;
+//    vTurbulenceNED(eY) = 0.0;
 
     break;
   }
-  case ttBerndt: {
+  case ttBerndt: { // This is very experimental and incomplete at the moment.
+
+    TurbGain = TurbGain * TurbGain * 100.0;
+  
     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));
@@ -411,7 +457,7 @@ void FGAtmosphere::Turbulence(void)
 
     vDirection.Normalize();
 
-    vTurbulence = TurbGain*Magnitude * vDirection;
+    vTurbulenceNED = TurbGain*Magnitude * vDirection;
     vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
 
     vBodyTurbGrad = Propagate->GetTl2b()*vTurbulenceGrad;
@@ -428,6 +474,56 @@ void FGAtmosphere::Turbulence(void)
 
     break;
   }
+  case ttCulp: { 
+
+    vTurbPQR(eP) = wind_from_clockwise;
+    if (TurbGain == 0.0) return;
+  
+    // keep the inputs within allowable limts for this model
+    if (TurbGain < 0.0) TurbGain = 0.0;
+    if (TurbGain > 1.0) TurbGain = 1.0;
+    if (TurbRate < 0.0) TurbRate = 0.0;
+    if (TurbRate > 30.0) TurbRate = 30.0;
+    if (Rhythmicity < 0.0) Rhythmicity = 0.0;
+    if (Rhythmicity > 1.0) Rhythmicity = 1.0;
+
+    // generate a sine wave corresponding to turbulence rate in hertz
+    double time = FDMExec->GetSimTime();
+    double sinewave = sin( time * TurbRate * 6.283185307 );
+
+    double random = 0.0;
+    if (target_time == 0.0) {
+      strength = random = 1 - 2.0*(double(rand())/double(RAND_MAX));
+      target_time = time + 0.71 + (random * 0.5);
+    }
+    if (time > target_time) {
+      spike = 1.0;
+      target_time = 0.0;
+    }    
+
+    // max vertical wind speed in fps, corresponds to TurbGain = 1.0
+    double max_vs = 40;
+
+    vTurbulenceNED(1) = vTurbulenceNED(2) = vTurbulenceNED(3) = 0.0;
+    double delta = strength * max_vs * TurbGain * (1-Rhythmicity) * spike;
+
+    // Vertical component of turbulence.
+    vTurbulenceNED(3) = sinewave * max_vs * TurbGain * Rhythmicity;
+    vTurbulenceNED(3)+= delta;
+    double HOverBMAC = Auxiliary->GetHOverBMAC();
+    if (HOverBMAC < 3.0)
+        vTurbulenceNED(3) *= HOverBMAC * 0.3333;
+    // Yaw component of turbulence.
+    vTurbulenceNED(1) = sin( delta * 3.0 );
+    vTurbulenceNED(2) = cos( delta * 3.0 );
+
+    // Roll component of turbulence. Clockwise vortex causes left roll.
+    vTurbPQR(eP) += delta * 0.04;
+
+    spike = spike * 0.9;
+    break;
+  }
   default:
     break;
   }
@@ -459,6 +555,7 @@ void FGAtmosphere::bind(void)
 {
   typedef double (FGAtmosphere::*PMF)(int) const;
   typedef double (FGAtmosphere::*PMFv)(void) const;
+  typedef void   (FGAtmosphere::*PMFd)(int,double);
   PropertyManager->Tie("atmosphere/T-R", this, (PMFv)&FGAtmosphere::GetTemperature);
   PropertyManager->Tie("atmosphere/rho-slugs_ft3", this, (PMFv)&FGAtmosphere::GetDensity);
   PropertyManager->Tie("atmosphere/P-psf", this, (PMFv)&FGAtmosphere::GetPressure);
@@ -475,34 +572,35 @@ void FGAtmosphere::bind(void)
   PropertyManager->Tie("atmosphere/delta-T", this, &FGAtmosphere::GetDeltaT, &FGAtmosphere::SetDeltaT);
   PropertyManager->Tie("atmosphere/T-sl-dev-F", this, &FGAtmosphere::GetSLTempDev, &FGAtmosphere::SetSLTempDev);
   PropertyManager->Tie("atmosphere/density-altitude", this, &FGAtmosphere::GetDensityAltitude);
+
+  PropertyManager->Tie("atmosphere/wind-north-fps", this, eNorth, (PMF)&FGAtmosphere::GetWindNED,
+                                                          (PMFd)&FGAtmosphere::SetWindNED);
+  PropertyManager->Tie("atmosphere/wind-east-fps",  this, eEast, (PMF)&FGAtmosphere::GetWindNED,
+                                                          (PMFd)&FGAtmosphere::SetWindNED);
+  PropertyManager->Tie("atmosphere/wind-down-fps",  this, eDown, (PMF)&FGAtmosphere::GetWindNED,
+                                                          (PMFd)&FGAtmosphere::SetWindNED);
+  PropertyManager->Tie("atmosphere/wind-from-cw", this, &FGAtmosphere::GetWindFromClockwise,
+                                                        &FGAtmosphere::SetWindFromClockwise);
+  PropertyManager->Tie("atmosphere/wind-mag-fps", this, &FGAtmosphere::GetWindspeed,
+                                                        &FGAtmosphere::SetWindspeed);
+  PropertyManager->Tie("atmosphere/total-wind-north-fps", this, eNorth, (PMF)&FGAtmosphere::GetTotalWindNED);
+  PropertyManager->Tie("atmosphere/total-wind-east-fps",  this, eEast,  (PMF)&FGAtmosphere::GetTotalWindNED);
+  PropertyManager->Tie("atmosphere/total-wind-down-fps",  this, eDown,  (PMF)&FGAtmosphere::GetTotalWindNED);
+
+  PropertyManager->Tie("atmosphere/gust-north-fps", this, eNorth, (PMF)&FGAtmosphere::GetGustNED,
+                                                          (PMFd)&FGAtmosphere::SetGustNED);
+  PropertyManager->Tie("atmosphere/gust-east-fps",  this, eEast, (PMF)&FGAtmosphere::GetGustNED,
+                                                          (PMFd)&FGAtmosphere::SetGustNED);
+  PropertyManager->Tie("atmosphere/gust-down-fps",  this, eDown, (PMF)&FGAtmosphere::GetGustNED,
+                                                          (PMFd)&FGAtmosphere::SetGustNED);
+
   PropertyManager->Tie("atmosphere/p-turb-rad_sec", this,1, (PMF)&FGAtmosphere::GetTurbPQR);
   PropertyManager->Tie("atmosphere/q-turb-rad_sec", this,2, (PMF)&FGAtmosphere::GetTurbPQR);
   PropertyManager->Tie("atmosphere/r-turb-rad_sec", this,3, (PMF)&FGAtmosphere::GetTurbPQR);
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-void FGAtmosphere::unbind(void)
-{
-  PropertyManager->Untie("atmosphere/T-R");
-  PropertyManager->Untie("atmosphere/rho-slugs_ft3");
-  PropertyManager->Untie("atmosphere/P-psf");
-  PropertyManager->Untie("atmosphere/a-fps");
-  PropertyManager->Untie("atmosphere/T-sl-R");
-  PropertyManager->Untie("atmosphere/rho-sl-slugs_ft3");
-  PropertyManager->Untie("atmosphere/P-sl-psf");
-  PropertyManager->Untie("atmosphere/a-sl-fps");
-  PropertyManager->Untie("atmosphere/delta-T");
-  PropertyManager->Untie("atmosphere/T-sl-dev-F");
-  PropertyManager->Untie("atmosphere/density-altitude");
-  PropertyManager->Untie("atmosphere/theta");
-  PropertyManager->Untie("atmosphere/sigma");
-  PropertyManager->Untie("atmosphere/delta");
-  PropertyManager->Untie("atmosphere/a-ratio");
-  PropertyManager->Untie("atmosphere/psiw-rad");
-  PropertyManager->Untie("atmosphere/p-turb-rad_sec");
-  PropertyManager->Untie("atmosphere/q-turb-rad_sec");
-  PropertyManager->Untie("atmosphere/r-turb-rad_sec");
+  PropertyManager->Tie("atmosphere/turb-rate", this, &FGAtmosphere::GetTurbRate, &FGAtmosphere::SetTurbRate);
+  PropertyManager->Tie("atmosphere/turb-gain", this, &FGAtmosphere::GetTurbGain, &FGAtmosphere::SetTurbGain);
+  PropertyManager->Tie("atmosphere/turb-rhythmicity", this, &FGAtmosphere::GetRhythmicity,
+                                                            &FGAtmosphere::SetRhythmicity);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -543,14 +641,16 @@ void FGAtmosphere::Debug(int from)
   if (debug_lvl & 16) { // Sanity checking
   }
   if (debug_lvl & 128) { // Turbulence
-    if (frame == 0 && from == 2) {
-      cout << "vTurbulence(X), vTurbulence(Y), vTurbulence(Z), "
+    if (first_pass && from == 2) {
+      first_pass = false;
+      cout << "vTurbulenceNED(X), vTurbulenceNED(Y), vTurbulenceNED(Z), "
            << "vTurbulenceGrad(X), vTurbulenceGrad(Y), vTurbulenceGrad(Z), "
            << "vDirection(X), vDirection(Y), vDirection(Z), "
            << "Magnitude, "
            << "vTurbPQR(P), vTurbPQR(Q), vTurbPQR(R), " << endl;
-    } else if (from == 2) {
-      cout << vTurbulence << ", " << vTurbulenceGrad << ", " << vDirection << ", " << Magnitude << ", " << vTurbPQR << endl;
+    } 
+    if (from == 2) {
+      cout << vTurbulenceNED << ", " << vTurbulenceGrad << ", " << vDirection << ", " << Magnitude << ", " << vTurbPQR << endl;
     }
   }
   if (debug_lvl & 64) {