]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/JSBSim/FGAtmosphere.cpp
Updated to match changes in radiostack.[ch]xx
[flightgear.git] / src / FDM / JSBSim / FGAtmosphere.cpp
index 085f95b694c9880bc67d270a8becc97062efaaae..e2ea8131d1aababaae5576e0a6f8ee43ea2ce97e 100644 (file)
@@ -36,6 +36,7 @@ HISTORY
 11/24/98   JSB   Created
 07/23/99   TP    Added implementation of 1959 Standard Atmosphere
                  Moved calculation of Mach number to FGTranslation
+                 Later updated to '76 model
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 COMMENTS, REFERENCES,  and NOTES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -56,10 +57,10 @@ INCLUDES
 #include "FGPosition.h"
 #include "FGAuxiliary.h"
 #include "FGOutput.h"
-#include "FGDefs.h"
 #include "FGMatrix33.h"
 #include "FGColumnVector3.h"
 #include "FGColumnVector4.h"
+#include "FGPropertyManager.h"
 
 static const char *IdSrc = "$Id$";
 static const char *IdHdr = ID_ATMOSPHERE;
@@ -69,12 +70,12 @@ CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 
-FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex),
-                                               vWindNED(3)
+FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
 {
   Name = "FGAtmosphere";
-  lastIndex=0;
-  h = 0;
+  lastIndex = 0;
+  h = 0.0;
+  psiw = 0.0;
   htab[0]=0;
   htab[1]=36089.239;
   htab[2]=65616.798;
@@ -84,25 +85,45 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex),
   htab[6]=200131.234;
   htab[7]=259186.352; //ft.
 
-  Calculate(h);
-  SLtemperature = temperature;
-  SLpressure    = pressure;
-  SLdensity     = density;
-  SLsoundspeed  = sqrt(SHRATIO*Reng*temperature);
-  rSLtemperature = 1.0/temperature;
-  rSLpressure    = 1.0/pressure;
-  rSLdensity     = 1.0/density;
-  rSLsoundspeed  = 1.0/SLsoundspeed;
-  useExternal=false;
-
-  if (debug_lvl & 2) cout << "Instantiated: " << Name << endl;
+  MagnitudedAccelDt = MagnitudeAccel = Magnitude = 0.0;
+  turbType = ttNone;
+//  turbType = ttBerndt; // temporarily disable turbulence until fully tested
+  TurbGain = 100.0;
+  
+  bind();
+  Debug(0);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 FGAtmosphere::~FGAtmosphere()
 {
-  if (debug_lvl & 2) cout << "Destroyed:    FGAtmosphere" << endl;
+  unbind();
+  Debug(1);
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+bool FGAtmosphere::InitModel(void)
+{
+  FGModel::InitModel();
+
+  Calculate(h);
+  SLtemperature = intTemperature;
+  SLpressure    = intPressure;
+  SLdensity     = intDensity;
+  SLsoundspeed  = sqrt(SHRatio*Reng*intTemperature);
+  rSLtemperature = 1.0/intTemperature;
+  rSLpressure    = 1.0/intPressure;
+  rSLdensity     = 1.0/intDensity;
+  rSLsoundspeed  = 1.0/SLsoundspeed;
+  temperature=&intTemperature;
+  pressure=&intPressure;
+  density=&intDensity;
+  
+  useExternal=false;
+  
+  return true;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -114,58 +135,58 @@ bool FGAtmosphere::Run(void)
     if (!useExternal) {
       h = Position->Geth();
       Calculate(h);
-    } else {
-      density = exDensity;
-      pressure = exPressure;
-      temperature = exTemperature;
+    } 
+
+    if (turbType != ttNone) {
+      Turbulence();
+      vWindNED += vTurbulence;
     }
 
     if (vWindNED(1) != 0.0) psiw = atan2( vWindNED(2), vWindNED(1) );
 
     if (psiw < 0) psiw += 2*M_PI;
 
-    soundspeed = sqrt(SHRATIO*Reng*temperature);
+    soundspeed = sqrt(SHRatio*Reng*(*temperature));
 
     State->Seta(soundspeed);
 
+    Debug(2);
+
   } else {                               // skip Run() execution this time
   }
+
   return false;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//
+// See reference 1
 
-void FGAtmosphere::Calculate(float altitude)
+void FGAtmosphere::Calculate(double altitude)
 {
-  //see reference [1]
+  double slope, reftemp, refpress;
+  int i = 0;
 
-  float slope,reftemp,refpress;
-  int i=0; bool lookup = false;
-  // cout << "Atmosphere:  h=" << altitude << " rho= " << density << endl;
-  i=lastIndex;
-  if(altitude < htab[lastIndex]) {
+  i = lastIndex;
+  if (altitude < htab[lastIndex]) {
     if (altitude <= 0) { 
-      i=0; altitude=0;
+      i = 0;
+      altitude=0;
     } else {
-       i=lastIndex-1;
-       while (htab[i] > altitude) { i--; }
+       i = lastIndex-1;
+       while (htab[i] > altitude) i--;
     }   
-  } else if (altitude > htab[lastIndex+1]){
-    if (altitude >= htab[7]){
-      i = 7; altitude = htab[7];
+  } else if (altitude > htab[lastIndex+1]) {
+    if (altitude >= htab[7]) {
+      i = 7;
+      altitude = htab[7];
     } else {
-      i=lastIndex+1;
-      while(htab[i+1] < altitude) { i++; } 
+      i = lastIndex+1;
+      while (htab[i+1] < altitude) i++;
     }  
   } 
 
   switch(i) {
-  case 0:     // sea level
-    slope     = -0.00356616; // R/ft.
-    reftemp   = 518.67;    // R
-    refpress  = 2116.22;    // psf
-    //refdens   = 0.00237767;  // slugs/cubic ft.
-    break;
   case 1:     // 36089 ft.
     slope     = 0;
     reftemp   = 389.97;
@@ -208,27 +229,205 @@ void FGAtmosphere::Calculate(float altitude)
     refpress  = 0.000122276;
     //refdens   = 2.19541e-10;
     break;
+  case 0:
+  default:     // sea level
+    slope     = -0.00356616; // R/ft.
+    reftemp   = 518.67;    // R
+    refpress  = 2116.22;    // psf
+    //refdens   = 0.00237767;  // slugs/cubic ft.
+    break;
+  
   }
  
   if (slope == 0) {
-    temperature = reftemp;
-    pressure = refpress*exp(-GRAVITY/(reftemp*Reng)*(altitude-htab[i]));
-    //density = refdens*exp(-GRAVITY/(reftemp*Reng)*(altitude-htab[i]));
-    density = pressure/(Reng*temperature);
+    intTemperature = reftemp;
+    intPressure = refpress*exp(-Inertial->SLgravity()/(reftemp*Reng)*(altitude-htab[i]));
+    //intDensity = refdens*exp(-Inertial->SLgravity()/(reftemp*Reng)*(altitude-htab[i]));
+    intDensity = intPressure/(Reng*intTemperature);
   } else {
-    temperature = reftemp+slope*(altitude-htab[i]);
-    pressure = refpress*pow(temperature/reftemp,-GRAVITY/(slope*Reng));
-    //density = refdens*pow(temperature/reftemp,-(GRAVITY/(slope*Reng)+1));
-    density = pressure/(Reng*temperature);
+    intTemperature = reftemp+slope*(altitude-htab[i]);
+    intPressure = refpress*pow(intTemperature/reftemp,-Inertial->SLgravity()/(slope*Reng));
+    //intDensity = refdens*pow(intTemperature/reftemp,-(Inertial->SLgravity()/(slope*Reng)+1));
+    intDensity = intPressure/(Reng*intTemperature);
   }
   lastIndex=i;
-  //cout << "Atmosphere:  h=" << altitude << " rho= " << density << endl;
+  //cout << "Atmosphere:  h=" << altitude << " rho= " << intDensity << endl;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-void FGAtmosphere::Debug(void)
+void FGAtmosphere::Turbulence(void)
 {
-    //TODO: Add your source code here
+  switch (turbType) {
+  case ttBerndt:
+    vDirectiondAccelDt(eX) = 1 - 2.0*(((double)(rand()))/RAND_MAX);
+    vDirectiondAccelDt(eY) = 1 - 2.0*(((double)(rand()))/RAND_MAX);
+    vDirectiondAccelDt(eZ) = 1 - 2.0*(((double)(rand()))/RAND_MAX);
+
+    MagnitudedAccelDt = 1 - 2.0*(((double)(rand()))/RAND_MAX);
+    MagnitudeAccel    += MagnitudedAccelDt*rate*State->Getdt();
+    Magnitude         += MagnitudeAccel*rate*State->Getdt();
+
+    vDirectiondAccelDt.Normalize();
+    vDirectionAccel += vDirectiondAccelDt*rate*State->Getdt();
+    vDirectionAccel.Normalize();
+    vDirection      += vDirectionAccel*rate*State->Getdt();
+    vDirection.Normalize();
+    
+    vTurbulence = TurbGain*Magnitude * vDirection;
+    vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
+
+    vBodyTurbGrad = State->GetTl2b()*vTurbulenceGrad;
+    vTurbPQR(eP) = vBodyTurbGrad(eY)/Aircraft->GetWingSpan();
+    if (Aircraft->GetHTailArm() != 0.0)
+      vTurbPQR(eQ) = vBodyTurbGrad(eZ)/Aircraft->GetHTailArm();
+    else
+      vTurbPQR(eQ) = vBodyTurbGrad(eZ)/10.0;
+
+    if (Aircraft->GetVTailArm())
+      vTurbPQR(eR) = vBodyTurbGrad(eX)/Aircraft->GetVTailArm();
+    else
+      vTurbPQR(eR) = vBodyTurbGrad(eX)/10.0;
+
+    break;
+  default:
+    break;
+  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void FGAtmosphere::UseExternal(void) {
+  temperature=&exTemperature;
+  pressure=&exPressure;
+  density=&exDensity;
+  useExternal=true;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void FGAtmosphere::UseInternal(void) {
+  temperature=&intTemperature;
+  pressure=&intPressure;
+  density=&intDensity;
+  useExternal=false;
+}
+  
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void FGAtmosphere::bind(void)
+{
+  typedef double (FGAtmosphere::*PMF)(int) const;
+  PropertyManager->Tie("atmosphere/T-R", this,
+                       &FGAtmosphere::GetTemperature);
+  PropertyManager->Tie("atmosphere/rho-slugs_ft3", this,
+                       &FGAtmosphere::GetDensity);
+  PropertyManager->Tie("atmosphere/P-psf", this,
+                       &FGAtmosphere::GetPressure);
+  PropertyManager->Tie("atmosphere/a-fps", this,
+                       &FGAtmosphere::GetSoundSpeed);
+  PropertyManager->Tie("atmosphere/T-sl-R", this,
+                       &FGAtmosphere::GetTemperatureSL);
+  PropertyManager->Tie("atmosphere/rho-sl-slugs_ft3", this,
+                       &FGAtmosphere::GetDensitySL);
+  PropertyManager->Tie("atmosphere/P-sl-psf", this,
+                       &FGAtmosphere::GetPressureSL);
+  PropertyManager->Tie("atmosphere/a-sl-fps", this,
+                       &FGAtmosphere::GetSoundSpeedSL);
+  PropertyManager->Tie("atmosphere/theta-norm", this,
+                       &FGAtmosphere::GetTemperatureRatio);
+  PropertyManager->Tie("atmosphere/sigma-norm", this,
+                       &FGAtmosphere::GetDensityRatio);
+  PropertyManager->Tie("atmosphere/delta-norm", this,
+                       &FGAtmosphere::GetPressureRatio);
+  PropertyManager->Tie("atmosphere/a-norm", this,
+                       &FGAtmosphere::GetSoundSpeedRatio);
+  PropertyManager->Tie("atmosphere/psiw-rad", this,
+                       &FGAtmosphere::GetWindPsi);
+  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/theta-norm");
+  PropertyManager->Untie("atmosphere/sigma-norm");
+  PropertyManager->Untie("atmosphere/delta-norm");
+  PropertyManager->Untie("atmosphere/a-norm");
+  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");
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//    The bitmasked value choices are as follows:
+//    unset: In this case (the default) JSBSim would only print
+//       out the normally expected messages, essentially echoing
+//       the config files as they are read. If the environment
+//       variable is not set, debug_lvl is set to 1 internally
+//    0: This requests JSBSim not to output any messages
+//       whatsoever.
+//    1: This value explicity requests the normal JSBSim
+//       startup messages
+//    2: This value asks for a message to be printed out when
+//       a class is instantiated
+//    4: When this value is set, a message is displayed when a
+//       FGModel object executes its Run() method
+//    8: When this value is set, various runtime state variables
+//       are printed out periodically
+//    16: When set various parameters are sanity checked and
+//       a message is printed out when they go out of bounds
+
+void FGAtmosphere::Debug(int from)
+{
+  if (debug_lvl <= 0) return;
+
+  if (debug_lvl & 1) { // Standard console startup message output
+    if (from == 0) { // Constructor
+    }
+  }
+  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
+    if (from == 0) cout << "Instantiated: FGAtmosphere" << endl;
+    if (from == 1) cout << "Destroyed:    FGAtmosphere" << endl;
+  }
+  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
+  }
+  if (debug_lvl & 8 ) { // Runtime state variables
+  }
+  if (debug_lvl & 16) { // Sanity checking
+  }
+  if (debug_lvl & 32) { // Turbulence
+    if (frame == 0 && from == 2) {
+      cout << "vTurbulence(X), vTurbulence(Y), vTurbulence(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 (debug_lvl & 64) {
+    if (from == 0) { // Constructor
+      cout << IdSrc << endl;
+      cout << IdHdr << endl;
+    }
+  }
 }