From 3cb45f498990d75b18e962af8f2b9bf09d09e45c Mon Sep 17 00:00:00 2001
From: ehofman <ehofman>
Date: Mon, 26 Jan 2009 20:36:17 +0000
Subject: [PATCH] Sync. with JSBSim CVS

---
 src/FDM/JSBSim/FGFDMExec.cpp                  |   3 +
 src/FDM/JSBSim/JSBSim.cxx                     |   4 +
 .../initialization/FGInitialCondition.cpp     |  13 +-
 src/FDM/JSBSim/input_output/FGScript.cpp      | 125 +++++++-------
 src/FDM/JSBSim/input_output/FGScript.h        |  16 +-
 src/FDM/JSBSim/math/FGColumnVector3.cpp       |   4 +-
 src/FDM/JSBSim/math/FGMatrix33.cpp            |   2 +-
 src/FDM/JSBSim/models/FGAircraft.cpp          |  17 +-
 src/FDM/JSBSim/models/FGAircraft.h            |   2 +-
 src/FDM/JSBSim/models/FGAtmosphere.cpp        |  78 +++++----
 src/FDM/JSBSim/models/FGAtmosphere.h          |   5 +
 src/FDM/JSBSim/models/FGAuxiliary.cpp         |   5 +
 src/FDM/JSBSim/models/FGAuxiliary.h           |  18 +-
 src/FDM/JSBSim/models/FGFCS.cpp               |   9 +-
 src/FDM/JSBSim/models/FGGroundReactions.cpp   |  19 +--
 src/FDM/JSBSim/models/FGInertial.cpp          |  14 +-
 src/FDM/JSBSim/models/FGLGear.cpp             | 155 +++++++++---------
 src/FDM/JSBSim/models/FGLGear.h               |   3 +
 src/FDM/JSBSim/models/FGOutput.cpp            |  28 ++--
 src/FDM/JSBSim/models/FGPropulsion.cpp        |   2 +-
 src/FDM/JSBSim/models/atmosphere/FGMSIS.cpp   |  12 +-
 src/FDM/JSBSim/models/propulsion/FGEngine.cpp |   3 +-
 src/FDM/JSBSim/models/propulsion/FGEngine.h   |   2 -
 src/FDM/JSBSim/models/propulsion/FGPiston.cpp |  86 +++-------
 src/FDM/JSBSim/models/propulsion/FGPiston.h   |   6 +-
 src/FDM/JSBSim/models/propulsion/FGRocket.cpp |   9 +-
 src/FDM/JSBSim/models/propulsion/FGTank.cpp   |   2 +
 .../JSBSim/models/propulsion/FGThruster.cpp   |   2 +-
 src/FDM/JSBSim/models/propulsion/FGThruster.h |   2 +-
 .../JSBSim/models/propulsion/FGTurbine.cpp    |  22 +--
 src/FDM/JSBSim/models/propulsion/FGTurbine.h  |   1 -
 31 files changed, 352 insertions(+), 317 deletions(-)

diff --git a/src/FDM/JSBSim/FGFDMExec.cpp b/src/FDM/JSBSim/FGFDMExec.cpp
index 2d5d409fa..25f849a9d 100644
--- a/src/FDM/JSBSim/FGFDMExec.cpp
+++ b/src/FDM/JSBSim/FGFDMExec.cpp
@@ -572,6 +572,9 @@ void FGFDMExec::BuildPropertyCatalog(struct PropertyCatalogStructure* pcs)
     sprintf(int_buf, "[%d]", node_idx);
     if (node_idx != 0) pcsNew->base_string += string(int_buf);
     if (pcs->node->getChild(i)->nChildren() == 0) {
+      if (pcsNew->base_string.substr(0,11) == string("/fdm/jsbsim")) {
+        pcsNew->base_string = pcsNew->base_string.erase(0,12);
+      }
       PropertyCatalog.push_back(pcsNew->base_string);
     } else {
       pcsNew->node = (FGPropertyManager*)pcs->node->getChild(i);
diff --git a/src/FDM/JSBSim/JSBSim.cxx b/src/FDM/JSBSim/JSBSim.cxx
index f63e106ad..a13821d82 100644
--- a/src/FDM/JSBSim/JSBSim.cxx
+++ b/src/FDM/JSBSim/JSBSim.cxx
@@ -282,6 +282,10 @@ FGJSBsim::FGJSBsim( double dt )
         fgGetDouble("/fdm/jsbsim/systems/hook/tailhook-offset-y-in", 0),
         fgGetDouble("/fdm/jsbsim/systems/hook/tailhook-offset-z-in", -16));
 
+    // Untie the write-state-file property to avoid creating an initfile.xml
+    // file on each reset.
+    fgGetNode("/fdm/jsbsim/simulation/write-state-file")->untie();
+
     crashed = false;
 }
 
diff --git a/src/FDM/JSBSim/initialization/FGInitialCondition.cpp b/src/FDM/JSBSim/initialization/FGInitialCondition.cpp
index eedec9345..2fbbae224 100644
--- a/src/FDM/JSBSim/initialization/FGInitialCondition.cpp
+++ b/src/FDM/JSBSim/initialization/FGInitialCondition.cpp
@@ -156,7 +156,13 @@ void FGInitialCondition::InitializeIC(void)
 
 void FGInitialCondition::WriteStateFile(int num)
 {
-  string filename = fdmex->GetFullAircraftPath() + "/" + "initfile.xml";
+  string filename = fdmex->GetFullAircraftPath();
+
+  if (filename.empty())
+    filename = "initfile.xml";
+  else
+    filename.append("/initfile.xml");
+  
   ofstream outfile(filename.c_str());
   FGPropagate* Propagate = fdmex->GetPropagate();
   
@@ -173,11 +179,10 @@ void FGInitialCondition::WriteStateFile(int num)
     outfile << "  <latitude unit=\"DEG\"> " << Propagate->GetLatitudeDeg() << " </latitude>" << endl;
     outfile << "  <altitude unit=\"FT\"> " << Propagate->Geth() << " </altitude>" << endl;
     outfile << "</initialize>" << endl;
+    outfile.close();
   } else {
-    cerr << "Could not open and/or write the state to the initial conditions file." << endl;
+    cerr << "Could not open and/or write the state to the initial conditions file: " << filename << endl;
   }
-
-  outfile.close();
 }
 
 //******************************************************************************
diff --git a/src/FDM/JSBSim/input_output/FGScript.cpp b/src/FDM/JSBSim/input_output/FGScript.cpp
index e0cddc34b..9ea8a9186 100755
--- a/src/FDM/JSBSim/input_output/FGScript.cpp
+++ b/src/FDM/JSBSim/input_output/FGScript.cpp
@@ -46,7 +46,6 @@ INCLUDES
 #include <initialization/FGTrim.h>
 
 #include <iostream>
-#include <iterator>
 
 namespace JSBSim {
 
@@ -204,12 +203,17 @@ bool FGScript::LoadScript( string script )
     // Retrieve the event name if given
     newEvent->Name = event_element->GetAttributeValue("name");
 
-    // Is this event persistent? That is, does it execute repeatedly as long as the
-    // condition is true, or does it execute as a one-shot event, only?
+    // Is this event persistent? That is, does it execute every time the
+    // condition triggers to true, or does it execute as a one-shot event, only?
     if (event_element->GetAttributeValue("persistent") == string("true")) {
       newEvent->Persistent = true;
     }
 
+    // Does this event execute continuously when triggered to true?
+    if (event_element->GetAttributeValue("continuous") == string("true")) {
+      newEvent->Continuous = true;
+    }
+
     // Process the conditions
     condition_element = event_element->FindElement("condition");
     if (condition_element != 0) {
@@ -304,7 +308,6 @@ bool FGScript::LoadScript( string script )
 
 bool FGScript::RunScript(void)
 {
-  vector <struct event>::iterator iEvent = Events.begin();
   unsigned i, j;
   unsigned event_ctr = 0;
 
@@ -314,90 +317,98 @@ bool FGScript::RunScript(void)
   if (currentTime > EndTime) return false; //Script done!
 
   // Iterate over all events.
-  while (iEvent < Events.end()) {
-    iEvent->PrevTriggered = iEvent->Triggered;
+  for (unsigned int ev_ctr=0; ev_ctr < Events.size(); ev_ctr++) {
     // Determine whether the set of conditional tests for this condition equate
-    // to true and should cause the event to execute.
-    if (iEvent->Condition->Evaluate()) {
-      if (!iEvent->Triggered) {
+    // to true and should cause the event to execute. If the conditions evaluate 
+    // to true, then the event is triggered. If the event is not persistent,
+    // then this trigger will remain set true. If the event is persistent,
+    // the trigger will reset to false when the condition evaluates to false.
+    if (Events[ev_ctr].Condition->Evaluate()) {
+      if (!Events[ev_ctr].Triggered) {
 
         // The conditions are true, do the setting of the desired Event parameters
-        for (i=0; i<iEvent->SetValue.size(); i++) {
-          iEvent->OriginalValue[i] = iEvent->SetParam[i]->getDoubleValue();
-          if (iEvent->Functions[i] != 0) { // Parameter should be set to a function value
-            iEvent->SetValue[i] = iEvent->Functions[i]->GetValue();
+        for (i=0; i<Events[ev_ctr].SetValue.size(); i++) {
+          Events[ev_ctr].OriginalValue[i] = Events[ev_ctr].SetParam[i]->getDoubleValue();
+          if (Events[ev_ctr].Functions[i] != 0) { // Parameter should be set to a function value
+            Events[ev_ctr].SetValue[i] = Events[ev_ctr].Functions[i]->GetValue();
           }
-          switch (iEvent->Type[i]) {
+          switch (Events[ev_ctr].Type[i]) {
           case FG_VALUE:
           case FG_BOOL:
-            iEvent->newValue[i] = iEvent->SetValue[i];
+            Events[ev_ctr].newValue[i] = Events[ev_ctr].SetValue[i];
             break;
           case FG_DELTA:
-            iEvent->newValue[i] = iEvent->OriginalValue[i] + iEvent->SetValue[i];
+            Events[ev_ctr].newValue[i] = Events[ev_ctr].OriginalValue[i] + Events[ev_ctr].SetValue[i];
             break;
           default:
             cerr << "Invalid Type specified" << endl;
             break;
           }
-          iEvent->StartTime = currentTime + iEvent->Delay;
-          iEvent->ValueSpan[i] = iEvent->newValue[i] - iEvent->OriginalValue[i];
-          iEvent->Transiting[i] = true;
+          Events[ev_ctr].StartTime = currentTime + Events[ev_ctr].Delay;
+          Events[ev_ctr].ValueSpan[i] = Events[ev_ctr].newValue[i] - Events[ev_ctr].OriginalValue[i];
+          Events[ev_ctr].Transiting[i] = true;
         }
       }
-      iEvent->Triggered = true;
-    } else if (iEvent->Persistent) {
-      iEvent->Triggered = false; // Reset the trigger for persistent events
-      iEvent->Notified = false;  // Also reset the notification flag
+      Events[ev_ctr].Triggered = true;
+
+    } else if (Events[ev_ctr].Persistent) { // If the event is persistent, reset the trigger.
+
+      Events[ev_ctr].Triggered = false; // Reset the trigger for persistent events
+      Events[ev_ctr].Notified = false;  // Also reset the notification flag
     }
 
-    if ((currentTime >= iEvent->StartTime) && iEvent->Triggered) {
+    if ((currentTime >= Events[ev_ctr].StartTime) && Events[ev_ctr].Triggered) {
 
-      for (i=0; i<iEvent->SetValue.size(); i++) {
-        if (iEvent->Transiting[i]) {
-          iEvent->TimeSpan = currentTime - iEvent->StartTime;
-          if (iEvent->Functions[i] == 0) {
-            switch (iEvent->Action[i]) {
-            case FG_RAMP:
-              if (iEvent->TimeSpan <= iEvent->TC[i]) {
-                newSetValue = iEvent->TimeSpan/iEvent->TC[i] * iEvent->ValueSpan[i] + iEvent->OriginalValue[i];
-              } else {
-                newSetValue = iEvent->newValue[i];
-                iEvent->Transiting[i] = false;
-              }
-              break;
-            case FG_STEP:
-              newSetValue = iEvent->newValue[i];
-              iEvent->Transiting[i] = false;
-              break;
-            case FG_EXP:
-              newSetValue = (1 - exp( -iEvent->TimeSpan/iEvent->TC[i] )) * iEvent->ValueSpan[i] + iEvent->OriginalValue[i];
-              break;
-            default:
-              cerr << "Invalid Action specified" << endl;
-              break;
+      for (i=0; i<Events[ev_ctr].SetValue.size(); i++) {
+        if (Events[ev_ctr].Transiting[i]) {
+          Events[ev_ctr].TimeSpan = currentTime - Events[ev_ctr].StartTime;
+          switch (Events[ev_ctr].Action[i]) {
+          case FG_RAMP:
+            if (Events[ev_ctr].TimeSpan <= Events[ev_ctr].TC[i]) {
+              newSetValue = Events[ev_ctr].TimeSpan/Events[ev_ctr].TC[i] * Events[ev_ctr].ValueSpan[i] + Events[ev_ctr].OriginalValue[i];
+            } else {
+              newSetValue = Events[ev_ctr].newValue[i];
+              if (Events[ev_ctr].Continuous != true) Events[ev_ctr].Transiting[i] = false;
             }
-          } else { // Set the new value based on a function
-            newSetValue = iEvent->Functions[i]->GetValue();
+            break;
+          case FG_STEP:
+            newSetValue = Events[ev_ctr].newValue[i];
+
+            // If this is not a continuous event, reset the transiting flag.
+            // Otherwise, it is known that the event is a continuous event.
+            // Furthermore, if the event is to be determined by a function,
+            // then the function will be continuously calculated.
+            if (Events[ev_ctr].Continuous != true)
+              Events[ev_ctr].Transiting[i] = false;
+            else if (Events[ev_ctr].Functions[i] != 0)
+              newSetValue = Events[ev_ctr].Functions[i]->GetValue();
+
+            break;
+          case FG_EXP:
+            newSetValue = (1 - exp( -Events[ev_ctr].TimeSpan/Events[ev_ctr].TC[i] )) * Events[ev_ctr].ValueSpan[i] + Events[ev_ctr].OriginalValue[i];
+            break;
+          default:
+            cerr << "Invalid Action specified" << endl;
+            break;
           }
-          iEvent->SetParam[i]->setDoubleValue(newSetValue);
+          Events[ev_ctr].SetParam[i]->setDoubleValue(newSetValue);
         }
       }
 
       // Print notification values after setting them
-      if (iEvent->Notify && !iEvent->Notified) {
-        cout << endl << "  Event " << event_ctr << " (" << iEvent->Name << ")"
+      if (Events[ev_ctr].Notify && !Events[ev_ctr].Notified) {
+        cout << endl << "  Event " << event_ctr << " (" << Events[ev_ctr].Name << ")"
              << " executed at time: " << currentTime << endl;
-        for (j=0; j<iEvent->NotifyProperties.size();j++) {
-          cout << "    " << iEvent->NotifyProperties[j]->GetName()
-               << " = " << iEvent->NotifyProperties[j]->getDoubleValue() << endl;
+        for (j=0; j<Events[ev_ctr].NotifyProperties.size();j++) {
+          cout << "    " << Events[ev_ctr].NotifyProperties[j]->GetName()
+               << " = " << Events[ev_ctr].NotifyProperties[j]->getDoubleValue() << endl;
         }
         cout << endl;
-        iEvent->Notified = true;
+        Events[ev_ctr].Notified = true;
       }
 
     }
 
-    iEvent++;
     event_ctr++;
   }
   return true;
diff --git a/src/FDM/JSBSim/input_output/FGScript.h b/src/FDM/JSBSim/input_output/FGScript.h
index f1a09bf1f..fdb0a97ff 100644
--- a/src/FDM/JSBSim/input_output/FGScript.h
+++ b/src/FDM/JSBSim/input_output/FGScript.h
@@ -70,13 +70,14 @@ CLASS DOCUMENTATION
     format. A test condition (or conditions) can be set up in an event in a
     script and when the condition evaluates to true, the specified
     action[s] is/are taken. An event can be <em>persistent</em>,
-    meaning that at all times when the test condition evaluates to true
-    the specified <em>set</em> actions take place. When the set of
-    tests evaluates to true for a given
+    meaning that at every time the test condition first evaluates to true
+    (toggling from false to true) then the specified <em>set</em> actions take
+    place. An event can also be defined to execute or evaluate continuously
+    while the condition is true. When the set of tests evaluates to true for a given
     condition, an item may be set to another value. This value may
     be a value, or a delta value, and the change from the
-    current value to the new value can be either via a step function,
-    a ramp, or an exponential approach. The speed of a ramp or
+    current value to the new value can be either via a step action,
+    a ramp, or an exponential approach. The speed of a ramp or exponential
     approach is specified via the time constant. Here is an example
     illustrating the format of the script file:
 
@@ -204,8 +205,8 @@ private:
   struct event {
     FGCondition     *Condition;
     bool             Persistent;
+    bool             Continuous;
     bool             Triggered;
-    bool             PrevTriggered;
     bool             Notify;
     bool             Notified;
     double           Delay;
@@ -226,8 +227,8 @@ private:
 
     event() {
       Triggered = false;
-      PrevTriggered = false;
       Persistent = false;
+      Continuous = false;
       Delay = 0.0;
       Notify = Notified = false;
       Name = "";
@@ -237,7 +238,6 @@ private:
 
     void reset(void) {
       Triggered = false;
-      PrevTriggered = false;
       Notified = false;
       StartTime = 0.0;
     }
diff --git a/src/FDM/JSBSim/math/FGColumnVector3.cpp b/src/FDM/JSBSim/math/FGColumnVector3.cpp
index 0bd5c011f..0bdeadaba 100644
--- a/src/FDM/JSBSim/math/FGColumnVector3.cpp
+++ b/src/FDM/JSBSim/math/FGColumnVector3.cpp
@@ -52,7 +52,7 @@ CLASS IMPLEMENTATION
 FGColumnVector3::FGColumnVector3(void)
 {
   data[0] = data[1] = data[2] = 0.0;
-  Debug(0);
+  // Debug(0);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -60,7 +60,7 @@ FGColumnVector3::FGColumnVector3(void)
 string FGColumnVector3::Dump(string delimeter) const
 {
   char buffer[256];
-  sprintf(buffer, "%10.3e%s%10.3e%s%10.3e", Entry(1), delimeter.c_str(), Entry(2), delimeter.c_str(), Entry(3));
+  sprintf(buffer, "%13.6e%s%13.6e%s%13.6e", Entry(1), delimeter.c_str(), Entry(2), delimeter.c_str(), Entry(3));
   return string(buffer);
 }
 
diff --git a/src/FDM/JSBSim/math/FGMatrix33.cpp b/src/FDM/JSBSim/math/FGMatrix33.cpp
index 9b4ce1231..91976f100 100644
--- a/src/FDM/JSBSim/math/FGMatrix33.cpp
+++ b/src/FDM/JSBSim/math/FGMatrix33.cpp
@@ -56,7 +56,7 @@ FGMatrix33::FGMatrix33(void)
   data[0] = data[1] = data[2] = data[3] = data[4] = data[5] =
     data[6] = data[7] = data[8] = 0.0;
 
-  Debug(0);
+  // Debug(0);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/src/FDM/JSBSim/models/FGAircraft.cpp b/src/FDM/JSBSim/models/FGAircraft.cpp
index 5dd8612c6..463e902a2 100644
--- a/src/FDM/JSBSim/models/FGAircraft.cpp
+++ b/src/FDM/JSBSim/models/FGAircraft.cpp
@@ -140,7 +140,7 @@ bool FGAircraft::Run(void)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-double FGAircraft::GetNlf(void)
+double FGAircraft::GetNlf(void) const
 {
   return -1*Aerodynamics->GetvFw(3)/MassBalance->GetWeight();
 }
@@ -216,13 +216,6 @@ void FGAircraft::bind(void)
   PropertyManager->Tie("metrics/lv-norm", this, &FGAircraft::Getlbarv);
   PropertyManager->Tie("metrics/vbarh-norm", this, &FGAircraft::Getvbarh);
   PropertyManager->Tie("metrics/vbarv-norm", this, &FGAircraft::Getvbarv);
-  PropertyManager->Tie("forces/hold-down", this, &FGAircraft::GetHoldDown, &FGAircraft::SetHoldDown);
-  PropertyManager->Tie("moments/l-total-lbsft", this, eL, (PMF)&FGAircraft::GetMoments);
-  PropertyManager->Tie("moments/m-total-lbsft", this, eM, (PMF)&FGAircraft::GetMoments);
-  PropertyManager->Tie("moments/n-total-lbsft", this, eN, (PMF)&FGAircraft::GetMoments);
-  PropertyManager->Tie("forces/fbx-total-lbs", this, eX, (PMF)&FGAircraft::GetForces);
-  PropertyManager->Tie("forces/fby-total-lbs", this, eY, (PMF)&FGAircraft::GetForces);
-  PropertyManager->Tie("forces/fbz-total-lbs", this, eZ, (PMF)&FGAircraft::GetForces);
   PropertyManager->Tie("metrics/aero-rp-x-in", this, eX, (PMF)&FGAircraft::GetXYZrp);
   PropertyManager->Tie("metrics/aero-rp-y-in", this, eY, (PMF)&FGAircraft::GetXYZrp);
   PropertyManager->Tie("metrics/aero-rp-z-in", this, eZ, (PMF)&FGAircraft::GetXYZrp);
@@ -232,6 +225,14 @@ void FGAircraft::bind(void)
   PropertyManager->Tie("metrics/visualrefpoint-x-in", this, eX, (PMF)&FGAircraft::GetXYZvrp);
   PropertyManager->Tie("metrics/visualrefpoint-y-in", this, eY, (PMF)&FGAircraft::GetXYZvrp);
   PropertyManager->Tie("metrics/visualrefpoint-z-in", this, eZ, (PMF)&FGAircraft::GetXYZvrp);
+  PropertyManager->Tie("forces/fbx-total-lbs", this, eX, (PMF)&FGAircraft::GetForces);
+  PropertyManager->Tie("forces/fby-total-lbs", this, eY, (PMF)&FGAircraft::GetForces);
+  PropertyManager->Tie("forces/fbz-total-lbs", this, eZ, (PMF)&FGAircraft::GetForces);
+  PropertyManager->Tie("forces/load-factor", this, &FGAircraft::GetNlf);
+  PropertyManager->Tie("forces/hold-down", this, &FGAircraft::GetHoldDown, &FGAircraft::SetHoldDown);
+  PropertyManager->Tie("moments/l-total-lbsft", this, eL, (PMF)&FGAircraft::GetMoments);
+  PropertyManager->Tie("moments/m-total-lbsft", this, eM, (PMF)&FGAircraft::GetMoments);
+  PropertyManager->Tie("moments/n-total-lbsft", this, eN, (PMF)&FGAircraft::GetMoments);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/src/FDM/JSBSim/models/FGAircraft.h b/src/FDM/JSBSim/models/FGAircraft.h
index f833c3a79..616da6cf0 100644
--- a/src/FDM/JSBSim/models/FGAircraft.h
+++ b/src/FDM/JSBSim/models/FGAircraft.h
@@ -171,7 +171,7 @@ public:
 
   void SetWingArea(double S) {WingArea = S;}
 
-  double GetNlf(void);
+  double GetNlf(void) const;
 
   inline FGColumnVector3& GetNwcg(void) { return vNwcg; }
 
diff --git a/src/FDM/JSBSim/models/FGAtmosphere.cpp b/src/FDM/JSBSim/models/FGAtmosphere.cpp
index 88ffbd5e0..215b3b8a1 100644
--- a/src/FDM/JSBSim/models/FGAtmosphere.cpp
+++ b/src/FDM/JSBSim/models/FGAtmosphere.cpp
@@ -71,13 +71,13 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
   h = 0.0;
   psiw = 0.0;
   htab[0]=0;
-  htab[1]=36089.239;
-  htab[2]=65616.798;
-  htab[3]=104986.878;
-  htab[4]=154199.475;
-  htab[5]=170603.675;
-  htab[6]=200131.234;
-  htab[7]=259186.352; //ft.
+  htab[1]= 36089.0;
+  htab[2]= 65617.0;
+  htab[3]=104987.0;
+  htab[4]=154199.0;
+  htab[5]=167322.0;
+  htab[6]=232940.0;
+  htab[7]=278385.0; //ft.
 
   MagnitudedAccelDt = MagnitudeAccel = Magnitude = 0.0;
   SetTurbType( ttCulp );
@@ -86,6 +86,8 @@ FGAtmosphere::FGAtmosphere(FGFDMExec* fdmex) : FGModel(fdmex)
   Rhythmicity = 0.1;
   spike = target_time = strength = 0.0;
   wind_from_clockwise = 0.0;
+  SutherlandConstant = 198.72; // deg Rankine
+  Beta = 2.269690E-08; // slug/(sec ft R^0.5)
 
   T_dev_sl = T_dev = delta_T = 0.0;
   StandardTempOnly = false;
@@ -174,52 +176,57 @@ void FGAtmosphere::Calculate(double altitude)
   }
 
   switch(i) {
-  case 1:     // 36089 ft.
+  case 0: // Sea level
+    slope     = -0.00356616; // R/ft.
+    reftemp   = 518.67;   // in degrees Rankine, 288.15 Kelvin
+    refpress  = 2116.22;    // psf
+    //refdens   = 0.00237767;  // slugs/cubic ft.
+    break;
+  case 1:     // 36089 ft. or 11 km
     slope     = 0;
-    reftemp   = 389.97;
-    refpress  = 472.452;
+    reftemp   = 389.97; // in degrees Rankine, 216.65 Kelvin
+    refpress  = 472.763;
     //refdens   = 0.000706032;
     break;
-  case 2:     // 65616 ft.
+  case 2:     // 65616 ft. or 20 km
     slope     = 0.00054864;
-    reftemp   = 389.97;
+    reftemp   = 389.97; // in degrees Rankine, 216.65 Kelvin
     refpress  = 114.636;
     //refdens   = 0.000171306;
     break;
-  case 3:     // 104986 ft.
-    slope     = 0.00153619;
-    reftemp   = 411.57;
-    refpress  = 8.36364;
+  case 3:     // 104986 ft. or 32 km
+    slope     = 0.001536192;
+    reftemp   = 411.57; // in degrees Rankine, 228.65 Kelvin
+    refpress  = 18.128;
     //refdens   = 1.18422e-05;
     break;
-  case 4:     // 154199 ft.
+  case 4:     // 154199 ft. 47 km
     slope     = 0;
-    reftemp   = 487.17;
-    refpress  = 0.334882;
+    reftemp   = 487.17; // in degrees Rankine, 270.65 Kelvin
+    refpress  = 2.316;
     //refdens   = 4.00585e-7;
     break;
-  case 5:     // 170603 ft.
-    slope     = -0.00109728;
-    reftemp   = 487.17;
-    refpress  = 0.683084;
+  case 5:     // 167322 ft. or 51 km
+    slope     = -0.001536192;
+    reftemp   = 487.17; // in degrees Rankine, 270.65 Kelvin
+    refpress  = 1.398;
     //refdens   = 8.17102e-7;
     break;
-  case 6:     // 200131 ft.
-    slope     = -0.00219456;
-    reftemp   = 454.17;
-    refpress  = 0.00684986;
+  case 6:     // 232940 ft. or 71 km
+    slope     = -0.00109728;
+    reftemp   = 386.368; // in degrees Rankine, 214.649 Kelvin
+    refpress  = 0.0826;
     //refdens   = 8.77702e-9;
     break;
-  case 7:     // 259186 ft.
+  case 7:     // 278385 ft. or 84.8520 km
     slope     = 0;
-    reftemp   = 325.17;
-    refpress  = 0.000122276;
+    reftemp   = 336.5; // in degrees Rankine, 186.94 Kelvin
+    refpress  = 0.00831;
     //refdens   = 2.19541e-10;
     break;
-  case 0:
   default:     // sea level
     slope     = -0.00356616; // R/ft.
-    reftemp   = 518.67;    // R
+    reftemp   = 518.67;   // in degrees Rankine, 288.15 Kelvin
     refpress  = 2116.22;    // psf
     //refdens   = 0.00237767;  // slugs/cubic ft.
     break;
@@ -250,7 +257,7 @@ void FGAtmosphere::Calculate(double altitude)
     intPressure = refpress*pow(intTemperature/reftemp,-Inertial->SLgravity()/(slope*Reng));
     intDensity = intPressure/(Reng*intTemperature);
   }
-
+  
   lastIndex=i;
 }
 
@@ -263,7 +270,7 @@ void FGAtmosphere::CalculateDerived(void)
   T_dev = (*temperature) - GetTemperature(h);
   density_altitude = h + T_dev * 66.7;
 
-  if (turbType == ttStandard || ttCulp) Turbulence();
+  if (turbType != ttNone) Turbulence();
 
   vTotalWindNED = vWindNED + vGustNED + vTurbulenceNED;
 
@@ -271,6 +278,9 @@ void FGAtmosphere::CalculateDerived(void)
   if (psiw < 0) psiw += 2*M_PI;
 
   soundspeed = sqrt(SHRatio*Reng*(*temperature));
+
+  intViscosity = Beta * pow(intTemperature, 1.5) / (SutherlandConstant + intTemperature);
+  intKinematicViscosity = intViscosity / intDensity;
 }
 
 
diff --git a/src/FDM/JSBSim/models/FGAtmosphere.h b/src/FDM/JSBSim/models/FGAtmosphere.h
index ed7bfcf93..b3f54f805 100644
--- a/src/FDM/JSBSim/models/FGAtmosphere.h
+++ b/src/FDM/JSBSim/models/FGAtmosphere.h
@@ -100,6 +100,10 @@ public:
   double GetDensity(double altitude);
   /// Returns the speed of sound in ft/sec.
   double GetSoundSpeed(void) const {return soundspeed;}
+  /// Returns the absolute viscosity.
+  double GetAbsoluteViscosity(void) const {return intViscosity;}
+  /// Returns the kinematic viscosity.
+  double GetKinematicViscosity(void) const {return intKinematicViscosity;}
 
   /// Returns the sea level temperature in degrees Rankine.
   double GetTemperatureSL(void) const { return SLtemperature; }
@@ -237,6 +241,7 @@ protected:
   bool useExternal;
   double exTemperature,exDensity,exPressure;
   double intTemperature, intDensity, intPressure;
+  double SutherlandConstant, Beta, intViscosity, intKinematicViscosity;
   double T_dev_sl, T_dev, delta_T, density_altitude;
   atmType atmosphere;
   bool StandardTempOnly;
diff --git a/src/FDM/JSBSim/models/FGAuxiliary.cpp b/src/FDM/JSBSim/models/FGAuxiliary.cpp
index a04ff725d..3fa74ead7 100755
--- a/src/FDM/JSBSim/models/FGAuxiliary.cpp
+++ b/src/FDM/JSBSim/models/FGAuxiliary.cpp
@@ -71,6 +71,7 @@ FGAuxiliary::FGAuxiliary(FGFDMExec* fdmex) : FGModel(fdmex)
   qbar = 0;
   qbarUW = 0.0;
   qbarUV = 0.0;
+  Re = 0.0;
   Mach = 0.0;
   alpha = beta = 0.0;
   adot = bdot = 0.0;
@@ -79,6 +80,7 @@ FGAuxiliary::FGAuxiliary(FGFDMExec* fdmex) : FGModel(fdmex)
   day_of_year = 1;
   seconds_in_day = 0.0;
   hoverbmac = hoverbcg = 0.0;
+  tatc = RankineToCelsius(tat);
 
   vPilotAccel.InitMatrix();
   vPilotAccelN.InitMatrix();
@@ -200,6 +202,8 @@ bool FGAuxiliary::Run()
     alpha = beta = adot = bdot = 0;
   }
 
+  Re = Vt * Aircraft->Getcbar() / Atmosphere->GetKinematicViscosity();
+
   qbar = 0.5*Atmosphere->GetDensity()*Vt*Vt;
   qbarUW = 0.5*Atmosphere->GetDensity()*(vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eW)*vAeroUVW(eW));
   qbarUV = 0.5*Atmosphere->GetDensity()*(vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eV)*vAeroUVW(eV));
@@ -349,6 +353,7 @@ void FGAuxiliary::bind(void)
   PropertyManager->Tie("aero/alpha-deg", this, inDegrees, (PMF)&FGAuxiliary::Getalpha);
   PropertyManager->Tie("aero/beta-deg", this, inDegrees, (PMF)&FGAuxiliary::Getbeta);
   PropertyManager->Tie("aero/mag-beta-deg", this, inDegrees, (PMF)&FGAuxiliary::GetMagBeta);
+  PropertyManager->Tie("aero/Re", this, &FGAuxiliary::GetReynoldsNumber);
   PropertyManager->Tie("aero/qbar-psf", this, &FGAuxiliary::Getqbar, &FGAuxiliary::Setqbar, true);
   PropertyManager->Tie("aero/qbarUW-psf", this, &FGAuxiliary::GetqbarUW, &FGAuxiliary::SetqbarUW, true);
   PropertyManager->Tie("aero/qbarUV-psf", this, &FGAuxiliary::GetqbarUV, &FGAuxiliary::SetqbarUV, true);
diff --git a/src/FDM/JSBSim/models/FGAuxiliary.h b/src/FDM/JSBSim/models/FGAuxiliary.h
index e551f5e3d..8d47ce8cc 100644
--- a/src/FDM/JSBSim/models/FGAuxiliary.h
+++ b/src/FDM/JSBSim/models/FGAuxiliary.h
@@ -181,14 +181,15 @@ public:
   double GetMagBeta (int unit) const { if (unit == inDegrees) return fabs(beta)*radtodeg;
                                        else cerr << "Bad units" << endl; return 0.0;}
 
-  double Getqbar    (void) const { return qbar;       }
-  double GetqbarUW  (void) const { return qbarUW;     }
-  double GetqbarUV  (void) const { return qbarUV;     }
-  double GetVt      (void) const { return Vt;         }
-  double GetVground (void) const { return Vground;    }
-  double GetMach    (void) const { return Mach;       }
-  double GetMachU   (void) const { return MachU;      }
-  double GetNz      (void) const { return Nz;         }
+  double Getqbar          (void) const { return qbar;       }
+  double GetqbarUW        (void) const { return qbarUW;     }
+  double GetqbarUV        (void) const { return qbarUV;     }
+  double GetReynoldsNumber(void) const { return Re;         }
+  double GetVt            (void) const { return Vt;         }
+  double GetVground       (void) const { return Vground;    }
+  double GetMach          (void) const { return Mach;       }
+  double GetMachU         (void) const { return MachU;      }
+  double GetNz            (void) const { return Nz;         }
 
   double GetHOverBCG(void) const { return hoverbcg; }
   double GetHOverBMAC(void) const { return hoverbmac; }
@@ -247,6 +248,7 @@ private:
 
   double Vt, Vground, Mach, MachU;
   double qbar, qbarUW, qbarUV;
+  double Re; // Reynolds Number = V*c/mu
   double alpha, beta;
   double adot,bdot;
   double psigt, gamma;
diff --git a/src/FDM/JSBSim/models/FGFCS.cpp b/src/FDM/JSBSim/models/FGFCS.cpp
index 95abe9ebb..3d974bf6e 100644
--- a/src/FDM/JSBSim/models/FGFCS.cpp
+++ b/src/FDM/JSBSim/models/FGFCS.cpp
@@ -46,7 +46,6 @@ INCLUDES
 #include <models/flight_control/FGDeadBand.h>
 #include <models/flight_control/FGGain.h>
 #include <models/flight_control/FGPID.h>
-#include <models/flight_control/FGGradient.h>
 #include <models/flight_control/FGSwitch.h>
 #include <models/flight_control/FGSummer.h>
 #include <models/flight_control/FGKinemat.h>
@@ -514,6 +513,10 @@ bool FGFCS::Load(Element* el, SystemType systype)
       return false;
     } else {
       document = LoadXMLDocument(file);
+      if (!document) {
+        cerr << "Error loading file " << file << endl;
+        return false;
+      }
       name = document->GetAttributeValue("name");
     }
   } else {
@@ -551,7 +554,7 @@ bool FGFCS::Load(Element* el, SystemType systype)
   }
 
   // After reading interface properties in a file, read properties in the local
-  // flight_control, autopiot, or system element. This allows general-purpose
+  // flight_control, autopilot, or system element. This allows general-purpose
   // systems to be defined in a file, with overrides or initial loaded constants
   // supplied in the relevant element of the aircraft configuration file.
 
@@ -967,7 +970,7 @@ void FGFCS::Debug(int from)
 
   if (debug_lvl & 1) { // Standard console startup message output
     if (from == 2) { // Loader
-      cout << endl << "  Flight Control (" << Name << ")" << endl;
+      cout << endl << "  " << Name << endl;
     }
   }
   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
diff --git a/src/FDM/JSBSim/models/FGGroundReactions.cpp b/src/FDM/JSBSim/models/FGGroundReactions.cpp
index f8067cfba..4638c1f39 100644
--- a/src/FDM/JSBSim/models/FGGroundReactions.cpp
+++ b/src/FDM/JSBSim/models/FGGroundReactions.cpp
@@ -89,17 +89,14 @@ bool FGGroundReactions::Run(void)
   vForces.InitMatrix();
   vMoments.InitMatrix();
 
-    // Sum forces and moments for all gear, here.
-    // Some optimizations may be made here - or rather in the gear code itself.
-    // The gear ::Run() method is called several times - once for each gear.
-    // Perhaps there is some commonality for things which only need to be
-    // calculated once.
-  if ( Propagate->GetDistanceAGL() < 300.0 ) { // Only execute gear code below 300 feet
-    for (unsigned int i=0; i<lGear.size(); i++) {
-      vForces  += lGear[i]->Force();
-      vMoments += lGear[i]->Moment();
-    }
-
+  // Sum forces and moments for all gear, here.
+  // Some optimizations may be made here - or rather in the gear code itself.
+  // The gear ::Run() method is called several times - once for each gear.
+  // Perhaps there is some commonality for things which only need to be
+  // calculated once.
+  for (unsigned int i=0; i<lGear.size(); i++) {
+    vForces  += lGear[i]->Force();
+    vMoments += lGear[i]->Moment();
   }
 
   return false;
diff --git a/src/FDM/JSBSim/models/FGInertial.cpp b/src/FDM/JSBSim/models/FGInertial.cpp
index 75a699f99..39aa91ddc 100644
--- a/src/FDM/JSBSim/models/FGInertial.cpp
+++ b/src/FDM/JSBSim/models/FGInertial.cpp
@@ -54,7 +54,7 @@ FGInertial::FGInertial(FGFDMExec* fgex) : FGModel(fgex)
 {
   Name = "FGInertial";
 
-  // Defaults
+  // Earth defaults
   RotationRate    = 0.00007292115;
   GM              = 14.07644180E15;     // WGS84 value
   RadiusReference = 20925650.00;        // Equatorial radius (WGS84)
@@ -64,6 +64,18 @@ FGInertial::FGInertial(FGFDMExec* fgex) : FGModel(fgex)
   b               = 20855486.5951;      // WGS84 semiminor axis length in feet
   earthPosAngle   = 0.0;
 
+  // Lunar defaults
+  /*
+  RotationRate    = 0.0000026617;
+  GM              = 1.7314079E14;         // Lunar GM
+  RadiusReference = 5702559.05;           // Equatorial radius
+  C2_0            = 0;                    // value for the C2,0 coefficient
+  J2              = 2.033542482111609E-4; // value for J2
+  a               = 5702559.05;           // semimajor axis length in feet 
+  b               = 5695439.63;           // semiminor axis length in feet
+  earthPosAngle   = 0.0;
+  */
+
   gAccelReference = GM/(RadiusReference*RadiusReference);
   gAccel          = GM/(RadiusReference*RadiusReference);
 
diff --git a/src/FDM/JSBSim/models/FGLGear.cpp b/src/FDM/JSBSim/models/FGLGear.cpp
index c1110dca6..080111976 100644
--- a/src/FDM/JSBSim/models/FGLGear.cpp
+++ b/src/FDM/JSBSim/models/FGLGear.cpp
@@ -57,8 +57,9 @@ static const char *IdHdr = ID_LGEAR;
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) : Exec(fdmex),
-                 GearNumber(number)
+FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
+  GearNumber(number),
+  Exec(fdmex)
 {
   Element *force_table=0;
   Element *dampCoeff=0;
@@ -341,97 +342,96 @@ FGColumnVector3& FGLGear::Force(void)
 
   if (isRetractable) ComputeRetractionState();
 
-  if (!GearDown) return vForce; // return the null vForce column vector
+  if (GearDown) {
 
-  vWhlBodyVec = MassBalance->StructuralToBody(vXYZ); // Get wheel in body frame
-  vLocalGear = Propagate->GetTb2l() * vWhlBodyVec; // Get local frame wheel location
+    vWhlBodyVec = MassBalance->StructuralToBody(vXYZ); // Get wheel in body frame
+    vLocalGear = Propagate->GetTb2l() * vWhlBodyVec; // Get local frame wheel location
 
-  gearLoc = Propagate->GetLocation().LocalToLocation(vLocalGear);
-  compressLength = -Exec->GetGroundCallback()->GetAGLevel(t, gearLoc, contact, normal, cvel);
+    gearLoc = Propagate->GetLocation().LocalToLocation(vLocalGear);
+    compressLength = -Exec->GetGroundCallback()->GetAGLevel(t, gearLoc, contact, normal, cvel);
 
-  // The compression length is measured in the Z-axis, only, at this time.
+    // The compression length is measured in the Z-axis, only, at this time.
 
-  if (compressLength > 0.00) {
+    if (compressLength > 0.00) {
 
-    WOW = true;
+      WOW = true;
 
-    // [The next equation should really use the vector to the contact patch of
-    // the tire including the strut compression and not the original vWhlBodyVec.]
+      // [The next equation should really use the vector to the contact patch of
+      // the tire including the strut compression and not the original vWhlBodyVec.]
 
-    vWhlVelVec      =  Propagate->GetTb2l() * (Propagate->GetPQR() * vWhlBodyVec);
-    vWhlVelVec     +=  Propagate->GetVel() - cvel;
-    compressSpeed   =  vWhlVelVec(eZ);
+      vWhlVelVec      =  Propagate->GetTb2l() * (Propagate->GetPQR() * vWhlBodyVec);
+      vWhlVelVec     +=  Propagate->GetVel() - cvel;
+      compressSpeed   =  vWhlVelVec(eZ);
 
-    InitializeReporting();
-    ComputeBrakeForceCoefficient();
-    ComputeSteeringAngle();
-    ComputeSlipAngle();
-    ComputeSideForceCoefficient();
-    ComputeVerticalStrutForce();
+      InitializeReporting();
+      ComputeBrakeForceCoefficient();
+      ComputeSteeringAngle();
+      ComputeSlipAngle();
+      ComputeSideForceCoefficient();
+      ComputeVerticalStrutForce();
 
-    // Compute the forces in the wheel ground plane.
+      // Compute the forces in the wheel ground plane.
 
-    double sign = RollingWhlVel>0?1.0:(RollingWhlVel<0?-1.0:0.0);
-    RollingForce = ((1.0 - TirePressureNorm) * 30 + vLocalForce(eZ) * BrakeFCoeff) * sign;
-    SideForce    = vLocalForce(eZ) * FCoeff;
+      double sign = RollingWhlVel>0?1.0:(RollingWhlVel<0?-1.0:0.0);
+      RollingForce = ((1.0 - TirePressureNorm) * 30 + vLocalForce(eZ) * BrakeFCoeff) * sign;
+      SideForce    = vLocalForce(eZ) * FCoeff;
 
-    // Transform these forces back to the local reference frame.
+      // Transform these forces back to the local reference frame.
 
-    vLocalForce(eX) = RollingForce*CosWheel - SideForce*SinWheel;
-    vLocalForce(eY) = SideForce*CosWheel    + RollingForce*SinWheel;
+      vLocalForce(eX) = RollingForce*CosWheel - SideForce*SinWheel;
+      vLocalForce(eY) = SideForce*CosWheel    + RollingForce*SinWheel;
 
-    // Transform the forces back to the body frame and compute the moment.
+      // Transform the forces back to the body frame and compute the moment.
 
-    vForce  = Propagate->GetTl2b() * vLocalForce;
+      vForce  = Propagate->GetTl2b() * vLocalForce;
 
-// Start experimental section for gear jitter reduction
-//
-// Lag and attenuate the XY-plane forces dependent on velocity
+      // Lag and attenuate the XY-plane forces dependent on velocity
 
-    double ca, cb, denom;
-    FGColumnVector3 Output;
+      double ca, cb, denom;
+      FGColumnVector3 Output;
 
-// This code implements a lag filter, C/(s + C) where
-// "C" is the filter coefficient. When "C" is chosen at the 
-// frame rate (in Hz), the jittering is significantly reduced. This is because
-// the jitter is present *at* the execution rate.
-// If a coefficient is set to something equal to or less than zero, the filter
-// is bypassed.
+      // This code implements a lag filter, C/(s + C) where
+      // "C" is the filter coefficient. When "C" is chosen at the 
+      // frame rate (in Hz), the jittering is significantly reduced. This is because
+      // the jitter is present *at* the execution rate.
+      // If a coefficient is set to something equal to or less than zero, the filter
+      // is bypassed.
 
-    if (LongForceLagFilterCoeff > 0) { 
-      denom = 2.00 + dT*LongForceLagFilterCoeff;
-      ca = dT*LongForceLagFilterCoeff / denom;
-      cb = (2.00 - dT*LongForceLagFilterCoeff) / denom;
-      Output(eX) = vForce(eX) * ca + prevIn(eX) * ca + prevOut(eX) * cb;
-      vForce(eX) = Output(eX);
+      if (LongForceLagFilterCoeff > 0) { 
+        denom = 2.00 + dT*LongForceLagFilterCoeff;
+        ca = dT*LongForceLagFilterCoeff / denom;
+        cb = (2.00 - dT*LongForceLagFilterCoeff) / denom;
+        Output(eX) = vForce(eX) * ca + prevIn(eX) * ca + prevOut(eX) * cb;
+        vForce(eX) = Output(eX);
+      }
+      if (LatForceLagFilterCoeff > 0) { 
+        denom = 2.00 + dT*LatForceLagFilterCoeff;
+        ca = dT*LatForceLagFilterCoeff / denom;
+        cb = (2.00 - dT*LatForceLagFilterCoeff) / denom;
+        Output(eY) = vForce(eY) * ca + prevIn(eY) * ca + prevOut(eY) * cb;
+        vForce(eY) = Output(eY);
+      }
+
+      prevIn = vForce;
+      prevOut = Output;
+
+      if ((fabs(RollingWhlVel) <= RFRV) && RFRV > 0) vForce(eX) *= fabs(RollingWhlVel)/RFRV;
+      if ((fabs(SideWhlVel) <= SFRV) && SFRV > 0) vForce(eY) *= fabs(SideWhlVel)/SFRV;
+
+  // End section for attentuating gear jitter
+
+      vMoment = vWhlBodyVec * vForce;
+
+    } else { // Gear is NOT compressed
+
+      WOW = false;
+      compressLength = 0.0;
+
+      // Return to neutral position between 1.0 and 0.8 gear pos.
+      SteerAngle *= max(GetGearUnitPos()-0.8, 0.0)/0.2;
+
+      ResetReporting();
     }
-    if (LatForceLagFilterCoeff > 0) { 
-      denom = 2.00 + dT*LatForceLagFilterCoeff;
-      ca = dT*LatForceLagFilterCoeff / denom;
-      cb = (2.00 - dT*LatForceLagFilterCoeff) / denom;
-      Output(eY) = vForce(eY) * ca + prevIn(eY) * ca + prevOut(eY) * cb;
-      vForce(eY) = Output(eY);
-    }
-
-    prevIn = vForce;
-    prevOut = Output;
-
-    if ((fabs(RollingWhlVel) <= RFRV) && RFRV > 0) vForce(eX) *= fabs(RollingWhlVel)/RFRV;
-    if ((fabs(SideWhlVel) <= SFRV) && SFRV > 0) vForce(eY) *= fabs(SideWhlVel)/SFRV;
-
-// End section for attentuating gear jitter
-
-    vMoment = vWhlBodyVec * vForce;
-
-  } else { // Gear is NOT compressed
-
-    WOW = false;
-    compressLength = 0.0;
-
-    // Return to neutral position between 1.0 and 0.8 gear pos.
-    SteerAngle *= max(GetGearUnitPos()-0.8, 0.0)/0.2;
-
-    ResetReporting();
   }
 
   ReportTakeoffOrLanding();
@@ -452,6 +452,7 @@ void FGLGear::ComputeRetractionState(void)
   double gearPos = GetGearUnitPos();
   if (gearPos < 0.01) {
     GearUp   = true;
+    WOW      = false;
     GearDown = false;
   } else if (gearPos > 0.99) {
     GearDown = true;
@@ -575,7 +576,7 @@ void FGLGear::ReportTakeoffOrLanding(void)
     LandingDistanceTraveled += Auxiliary->GetVground()*deltaT;
 
   if (StartedGroundRun) {
-     TakeoffDistanceTraveled50ft += Auxiliary->GetVground()*deltaT;
+    TakeoffDistanceTraveled50ft += Auxiliary->GetVground()*deltaT;
     if (WOW) TakeoffDistanceTraveled += Auxiliary->GetVground()*deltaT;
   }
 
@@ -739,6 +740,8 @@ void FGLGear::bind(void)
     snprintf(property_name, 80, "gear/unit[%d]/z-position", GearNumber);
     Exec->GetPropertyManager()->Tie( property_name, (FGLGear*)this,
                           &FGLGear::GetZPosition, &FGLGear::SetZPosition);
+    snprintf(property_name, 80, "gear/unit[%d]/compression-ft", GearNumber);
+    Exec->GetPropertyManager()->Tie( property_name, &compressLength );
   }
 
   if( isRetractable ) {
@@ -752,6 +755,8 @@ void FGLGear::bind(void)
 
 void FGLGear::Report(ReportType repType)
 {
+  if (fabs(TakeoffDistanceTraveled) < 0.001) return; // Don't print superfluous reports
+
   switch(repType) {
   case erLand:
     cout << endl << "Touchdown report for " << name << endl;
diff --git a/src/FDM/JSBSim/models/FGLGear.h b/src/FDM/JSBSim/models/FGLGear.h
index 3dc81b7c1..4de980edf 100644
--- a/src/FDM/JSBSim/models/FGLGear.h
+++ b/src/FDM/JSBSim/models/FGLGear.h
@@ -257,6 +257,9 @@ public:
   /// Sets the brake value in percent (0 - 100)
   inline void SetBrake(double bp) {brakePct = bp;}
 
+  /// Sets the weight-on-wheels flag.
+  void SetWOW(bool wow) {WOW = wow;}
+
   /** Set the console touchdown reporting feature
       @param flag true turns on touchdown reporting, false turns it off */
   inline void SetReport(bool flag) { ReportEnable = flag; }
diff --git a/src/FDM/JSBSim/models/FGOutput.cpp b/src/FDM/JSBSim/models/FGOutput.cpp
index 032245e85..4b3a95f88 100644
--- a/src/FDM/JSBSim/models/FGOutput.cpp
+++ b/src/FDM/JSBSim/models/FGOutput.cpp
@@ -247,6 +247,7 @@ void FGOutput::DelimitedOutput(string fname)
     if (SubSystems & ssVelocities) {
       outstream << delimeter;
       outstream << "q bar (psf)" + delimeter;
+      outstream << "Reynolds Number" + delimeter;
       outstream << "V_{Total} (ft/s)" + delimeter;
       outstream << "V_{Inertial} (ft/s)" + delimeter;
       outstream << "UBody" + delimeter + "VBody" + delimeter + "WBody" + delimeter;
@@ -266,6 +267,9 @@ void FGOutput::DelimitedOutput(string fname)
     if (SubSystems & ssAtmosphere) {
       outstream << delimeter;
       outstream << "Rho (slugs/ft^3)" + delimeter;
+      outstream << "Absolute Viscosity" + delimeter;
+      outstream << "Kinematic Viscosity" + delimeter;
+      outstream << "Temperature (R)" + delimeter;
       outstream << "P_{SL} (psf)" + delimeter;
       outstream << "P_{Ambient} (psf)" + delimeter;
       outstream << "Turbulence Magnitude (ft/sec)" + delimeter;
@@ -274,17 +278,17 @@ void FGOutput::DelimitedOutput(string fname)
     }
     if (SubSystems & ssMassProps) {
       outstream << delimeter;
-      outstream << "I_xx" + delimeter;
-      outstream << "I_xy" + delimeter;
-      outstream << "I_xz" + delimeter;
-      outstream << "I_yx" + delimeter;
-      outstream << "I_yy" + delimeter;
-      outstream << "I_yz" + delimeter;
-      outstream << "I_zx" + delimeter;
-      outstream << "I_zy" + delimeter;
-      outstream << "I_zz" + delimeter;
+      outstream << "I_{xx}" + delimeter;
+      outstream << "I_{xy}" + delimeter;
+      outstream << "I_{xz}" + delimeter;
+      outstream << "I_{yx}" + delimeter;
+      outstream << "I_{yy}" + delimeter;
+      outstream << "I_{yz}" + delimeter;
+      outstream << "I_{zx}" + delimeter;
+      outstream << "I_{zy}" + delimeter;
+      outstream << "I_{zz}" + delimeter;
       outstream << "Mass" + delimeter;
-      outstream << "X_cg" + delimeter + "Y_cg" + delimeter + "Z_cg";
+      outstream << "X_{cg}" + delimeter + "Y_{cg}" + delimeter + "Z_{cg}";
     }
     if (SubSystems & ssPropagate) {
       outstream << delimeter;
@@ -348,6 +352,7 @@ void FGOutput::DelimitedOutput(string fname)
   if (SubSystems & ssVelocities) {
     outstream << delimeter;
     outstream << Auxiliary->Getqbar() << delimeter;
+    outstream << Auxiliary->GetReynoldsNumber() << delimeter;
     outstream << setprecision(12) << Auxiliary->GetVt() << delimeter;
     outstream << Propagate->GetInertialVelocityMagnitude() << delimeter;
     outstream << setprecision(12) << Propagate->GetUVW().Dump(delimeter) << delimeter;
@@ -367,6 +372,9 @@ void FGOutput::DelimitedOutput(string fname)
   if (SubSystems & ssAtmosphere) {
     outstream << delimeter;
     outstream << Atmosphere->GetDensity() << delimeter;
+    outstream << Atmosphere->GetAbsoluteViscosity() << delimeter;
+    outstream << Atmosphere->GetKinematicViscosity() << delimeter;
+    outstream << Atmosphere->GetTemperature() << delimeter;
     outstream << Atmosphere->GetPressureSL() << delimeter;
     outstream << Atmosphere->GetPressure() << delimeter;
     outstream << Atmosphere->GetTurbMagnitude() << delimeter;
diff --git a/src/FDM/JSBSim/models/FGPropulsion.cpp b/src/FDM/JSBSim/models/FGPropulsion.cpp
index 2ec2e2bb3..82092f8ad 100644
--- a/src/FDM/JSBSim/models/FGPropulsion.cpp
+++ b/src/FDM/JSBSim/models/FGPropulsion.cpp
@@ -184,7 +184,7 @@ bool FGPropulsion::GetSteadyState(void)
       while (!steady && j < 6000) {
         Engines[i]->Calculate();
         lastThrust = currentThrust;
-        currentThrust = Engines[i]->GetThrust();
+        currentThrust = Engines[i]->GetThruster()->GetThrust();
         if (fabs(lastThrust-currentThrust) < 0.0001) {
           steady_count++;
           if (steady_count > 120) {
diff --git a/src/FDM/JSBSim/models/atmosphere/FGMSIS.cpp b/src/FDM/JSBSim/models/atmosphere/FGMSIS.cpp
index 6eeff5af8..e3197404d 100755
--- a/src/FDM/JSBSim/models/atmosphere/FGMSIS.cpp
+++ b/src/FDM/JSBSim/models/atmosphere/FGMSIS.cpp
@@ -136,7 +136,7 @@ bool MSIS::InitModel(void)
   pressure = &intPressure;
   density = &intDensity;
 
-  useExternal=false;
+  UseInternal();
 
   return true;
 }
@@ -174,19 +174,11 @@ bool MSIS::Run(void)
     intTemperature = output.t[1] * 1.8;
     intDensity     = output.d[5] * 1.940321;
     intPressure    = 1716.488 * intDensity * intTemperature;
-    soundspeed     = sqrt(2403.0832 * intTemperature);
     //cout << "T=" << intTemperature << " D=" << intDensity << " P=";
     //cout << intPressure << " a=" << soundspeed << endl;
   }
 
-  if (turbType != ttNone) {
-    Turbulence();
-    vWindNED += vTurbulenceNED;
-  }
-
-  if (vWindNED(1) != 0.0) psiw = atan2( vWindNED(2), vWindNED(1) );
-
-  if (psiw < 0) psiw += 2*M_PI;
+  CalculateDerived();
 
   Debug(2);
 
diff --git a/src/FDM/JSBSim/models/propulsion/FGEngine.cpp b/src/FDM/JSBSim/models/propulsion/FGEngine.cpp
index 49f67e2f4..49937ab8f 100644
--- a/src/FDM/JSBSim/models/propulsion/FGEngine.cpp
+++ b/src/FDM/JSBSim/models/propulsion/FGEngine.cpp
@@ -119,7 +119,7 @@ FGEngine::FGEngine(FGFDMExec* exec, Element* engine_element, int engine_number)
   snprintf(property_name, 80, "propulsion/engine[%d]/set-running", EngineNumber);
   PropertyManager->Tie( property_name, this, &FGEngine::GetRunning, &FGEngine::SetRunning );
   snprintf(property_name, 80, "propulsion/engine[%u]/thrust-lbs", EngineNumber);
-  PropertyManager->Tie( property_name, this, &FGEngine::GetThrust);
+  PropertyManager->Tie( property_name, Thruster, &FGThruster::GetThrust);
   snprintf(property_name, 80, "propulsion/engine[%u]/fuel-flow-rate-pps", EngineNumber);
   PropertyManager->Tie( property_name, this, &FGEngine::GetFuelFlowRate);
 
@@ -138,7 +138,6 @@ FGEngine::~FGEngine()
 
 void FGEngine::ResetToIC(void)
 {
-  Thrust = 0.0;
   Throttle = 0.0;
   Mixture = 1.0;
   Starter = false;
diff --git a/src/FDM/JSBSim/models/propulsion/FGEngine.h b/src/FDM/JSBSim/models/propulsion/FGEngine.h
index dd8b1302b..783fbc0fe 100644
--- a/src/FDM/JSBSim/models/propulsion/FGEngine.h
+++ b/src/FDM/JSBSim/models/propulsion/FGEngine.h
@@ -150,7 +150,6 @@ public:
   virtual double getFuelFlow_gph () const {return FuelFlow_gph;}
   virtual double getFuelFlow_pph () const {return FuelFlow_pph;}
   virtual double GetFuelFlowRate(void) const {return FuelFlowRate;}
-  virtual double GetThrust(void) const { return Thrust; }
   virtual bool   GetStarved(void) { return Starved; }
   virtual bool   GetRunning(void) const { return Running; }
   virtual bool   GetCranking(void) { return Cranking; }
@@ -216,7 +215,6 @@ protected:
   double MaxThrottle;
   double MinThrottle;
 
-  double Thrust;
   double Throttle;
   double Mixture;
   double FuelExpended;
diff --git a/src/FDM/JSBSim/models/propulsion/FGPiston.cpp b/src/FDM/JSBSim/models/propulsion/FGPiston.cpp
index e23d981a8..bac4eaa91 100644
--- a/src/FDM/JSBSim/models/propulsion/FGPiston.cpp
+++ b/src/FDM/JSBSim/models/propulsion/FGPiston.cpp
@@ -132,21 +132,6 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number)
   *Lookup_Combustion_Efficiency << 1.60 << 0.525;
   *Lookup_Combustion_Efficiency << 2.00 << 0.345;
 
-  Power_Mixture_Correlation = new FGTable(13);
-  *Power_Mixture_Correlation << (14.7/1.6) << 0.780;
-  *Power_Mixture_Correlation << 10 <<  0.860;
-  *Power_Mixture_Correlation << 11 <<  0.935;
-  *Power_Mixture_Correlation << 12 <<  0.980;
-  *Power_Mixture_Correlation << 13 <<  1.000;
-  *Power_Mixture_Correlation << 14 <<  0.990;
-  *Power_Mixture_Correlation << 15 <<  0.964;
-  *Power_Mixture_Correlation << 16 <<  0.925;
-  *Power_Mixture_Correlation << 17 <<  0.880;
-  *Power_Mixture_Correlation << 18 <<  0.830;
-  *Power_Mixture_Correlation << 19 <<  0.785;
-  *Power_Mixture_Correlation << 20 <<  0.740;
-  *Power_Mixture_Correlation << (14.7/0.6) << 0.58;
-
   Mixture_Efficiency_Correlation = new FGTable(15);
   *Mixture_Efficiency_Correlation << 0.05000 << 0.00000;
   *Mixture_Efficiency_Correlation << 0.05137 << 0.00862;
@@ -165,21 +150,6 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number)
   *Mixture_Efficiency_Correlation << 0.12500 << 0.00000;
 
 
-/*
-Manifold_Pressure_Lookup = new
-
-        0       0.2      0.4      0.6      0.8     1
-0        1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
-1000    0.7778    0.8212    0.8647    0.9081    0.9516    0.9950
-2000    0.5556    0.6424    0.7293    0.8162    0.9031    0.9900
-3000    0.3333    0.4637    0.5940    0.7243    0.8547    0.9850
-4000    0.2000    0.2849    0.4587    0.6324    0.8062    0.9800
-5000    0.2000    0.2000    0.3233    0.5406    0.7578    0.9750
-6000    0.2000    0.2000    0.2000    0.4487    0.7093    0.9700
-7000    0.2000    0.2000    0.2000    0.2000    0.4570    0.7611
-8000    0.2000    0.2000    0.2000    0.2000    0.2047    0.5522
-*/
-
   // Read inputs from engine data file where present.
 
   if (el->FindElement("minmp")) // Should have ELSE statement telling default value used?
@@ -238,10 +208,17 @@ Manifold_Pressure_Lookup = new
       RatedAltitude[2] = el->FindElementValueAsNumberConvertTo("ratedaltitude3", "FT");
   }
 
+  MaxManifoldPressure_Percent = MaxManifoldPressure_inHg / 29.92;
   // Create a BSFC to match the engine if not provided
   if (BSFC < 0) {
       BSFC = ( Displacement * MaxRPM * volumetric_efficiency ) / (9411 * MaxHP);
+      BSFC *= (MaxManifoldPressure_Percent * MaxManifoldPressure_Percent * MaxManifoldPressure_Percent);
   }
+  if ( MaxManifoldPressure_inHg > 29.9 ) {   // Don't allow boosting with a bogus number
+      MaxManifoldPressure_inHg = 29.9;
+      MaxManifoldPressure_Percent = MaxManifoldPressure_inHg / 29.92;
+  }
+
   char property_name[80];
   snprintf(property_name, 80, "propulsion/engine[%d]/power-hp", EngineNumber);
   PropertyManager->Tie(property_name, &HP);
@@ -299,7 +276,9 @@ Manifold_Pressure_Lookup = new
     BoostSpeed = 0;
   }
   bBoostOverride = (BoostOverride == 1 ? true : false);
-  if (MinThrottle < 0.001) MinThrottle = 0.001;  //MinThrottle is a denominator in a power equation so it can't be zero
+  if (MinThrottle < 0.12) MinThrottle = 0.12;  //MinThrottle is limited to 0.12 to prevent the
+                                               // throttle area equation from going negative
+                                               // 0.12 is 1% of maximum area
   Debug(0); // Call Debug() routine from constructor if needed
 }
 
@@ -308,7 +287,6 @@ Manifold_Pressure_Lookup = new
 FGPiston::~FGPiston()
 {
   delete Lookup_Combustion_Efficiency;
-  delete Power_Mixture_Correlation;
   delete Mixture_Efficiency_Correlation;
   Debug(1); // Call Debug() routine from constructor if needed
 }
@@ -338,7 +316,8 @@ double FGPiston::Calculate(void)
   if (FuelFlow_gph > 0.0) ConsumeFuel();
 
   Throttle = FCS->GetThrottlePos(EngineNumber);
-  ThrottlePos = MinThrottle+((MaxThrottle-MinThrottle)*Throttle );
+  // calculate the throttle plate angle.  1 unit is pi/2 radians.
+  ThrottleAngle = MinThrottle+((MaxThrottle-MinThrottle)*Throttle );
   Mixture = FCS->GetMixturePos(EngineNumber);
 
   //
@@ -367,8 +346,7 @@ double FGPiston::Calculate(void)
 //    Running = false;
 
   doEnginePower();
-if(HP<0.1250)
-  Running = false;
+  if (HP < 0.1250) Running = false;
 
   doEGT();
   doCHT();
@@ -390,8 +368,6 @@ if(HP<0.1250)
 double FGPiston::CalcFuelNeed(void)
 {
   double dT = State->Getdt() * Propulsion->GetRate();
-  FuelFlow_pph = FuelFlow_gph * 6.0; // Assumes 6 lbs / gallon
-  FuelFlowRate = FuelFlow_pph / 3600.0;
   FuelExpended = FuelFlowRate * dT;
   return FuelExpended;
 }
@@ -510,20 +486,18 @@ void FGPiston::doBoostControl(void)
  * from the throttle position, turbo/supercharger boost control
  * system, engine speed and local ambient air density.
  *
- * TODO: changes in MP should not be instantaneous -- introduce
- * a lag between throttle changes and MP changes, to allow pressure
- * to build up or disperse.
- *
- * Inputs: minMAP, maxMAP, p_amb, Throttle
+ * Inputs: p_amb, Throttle, MaxManifoldPressure_Percent, ThrottleAngle
+ *         RPM, MaxRPM
  *
  * Outputs: MAP, ManifoldPressure_inHg
  */
 
 void FGPiston::doMAP(void)
 {
-    suction_loss = RPM > 0.0 ? ThrottlePos * MaxRPM / RPM : 1.0;
-    if (suction_loss > 1.0) suction_loss = 1.0;
-    MAP = p_amb * suction_loss;
+ // estimate throttle plate area.  This maps 0.2 -> 0.1 for historical performance reasons
+    double throttle_area = ThrottleAngle * 1.125 - 0.125;
+    map_coefficient = pow ((throttle_area * MaxManifoldPressure_Percent),RPM/MaxRPM);
+    MAP = p_amb * map_coefficient;
 
     if(Boosted) {
       // If takeoff boost is fitted, we currently assume the following throttle map:
@@ -549,7 +523,7 @@ void FGPiston::doMAP(void)
         }
       }
       // Boost the manifold pressure.
-      double boost_factor = BoostMul[BoostSpeed] * suction_loss * RPM/RatedRPM[BoostSpeed];
+      double boost_factor = BoostMul[BoostSpeed] * map_coefficient * RPM/RatedRPM[BoostSpeed];
       if (boost_factor < 1.0) boost_factor = 1.0;  // boost will never reduce the MAP
       MAP *= boost_factor;
       // Now clip the manifold pressure to BCV or Wastegate setting.
@@ -575,7 +549,7 @@ void FGPiston::doMAP(void)
  * (used in CHT calculation for air-cooled engines).
  *
  * Inputs: p_amb, R_air, T_amb, MAP, Displacement,
- *   RPM, volumetric_efficiency, ThrottlePos
+ *   RPM, volumetric_efficiency, ThrottleAngle
  *
  * TODO: Model inlet manifold air temperature.
  *
@@ -587,7 +561,7 @@ void FGPiston::doAirFlow(void)
   rho_air = p_amb / (R_air * T_amb);
   double displacement_SI = Displacement * in3tom3;
   double swept_volume = (displacement_SI * (RPM/60)) / 2;
-  double v_dot_air = swept_volume * volumetric_efficiency * suction_loss;
+  double v_dot_air = swept_volume * volumetric_efficiency * map_coefficient;
 
   double rho_air_manifold = MAP / (R_air * T_amb);
   m_dot_air = v_dot_air * rho_air_manifold;
@@ -609,11 +583,9 @@ void FGPiston::doFuelFlow(void)
 //  double AFR = 10+(12*(1-Mixture));// mixture 10:1 to 22:1
 //  m_dot_fuel = m_dot_air / AFR;
   m_dot_fuel = (m_dot_air * equivalence_ratio) / 14.7;
-  FuelFlow_gph = m_dot_fuel
-    * 3600            // seconds to hours
-    * 2.2046            // kg to lb
-    / 6.0;            // lb to gal_us of gasoline
-//    / 6.6;            // lb to gal_us of kerosene
+  FuelFlowRate =  m_dot_fuel * 2.2046;  // kg to lb
+  FuelFlow_pph = FuelFlowRate  * 3600;  // seconds to hours
+  FuelFlow_gph = FuelFlow_pph / 6.0;    // Assumes 6 lbs / gallon
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -649,7 +621,7 @@ void FGPiston::doEnginePower(void)
     if ( Magnetos != 3 ) power *= SparkFailDrop;
 
 
-    HP = (FuelFlow_gph * 6.0 / BSFC )* ME * suction_loss * power;
+    HP = (FuelFlow_gph * 6.0 / BSFC )* ME * map_coefficient * power;
 
   } else {
 
@@ -871,17 +843,13 @@ void FGPiston::Debug(int from)
       cout << "      IdleRPM: "             << IdleRPM                  << endl;
       cout << "      MaxThrottle: "         << MaxThrottle              << endl;
       cout << "      MinThrottle: "         << MinThrottle              << endl;
+      cout << "      BSFC: "                << BSFC                     << endl;
 
       cout << endl;
       cout << "      Combustion Efficiency table:" << endl;
       Lookup_Combustion_Efficiency->Print();
       cout << endl;
 
-      cout << endl;
-      cout << "      Power Mixture Correlation table:" << endl;
-      Power_Mixture_Correlation->Print();
-      cout << endl;
-
       cout << endl;
       cout << "      Mixture Efficiency Correlation table:" << endl;
       Mixture_Efficiency_Correlation->Print();
diff --git a/src/FDM/JSBSim/models/propulsion/FGPiston.h b/src/FDM/JSBSim/models/propulsion/FGPiston.h
index d4fa080f3..8878c5040 100644
--- a/src/FDM/JSBSim/models/propulsion/FGPiston.h
+++ b/src/FDM/JSBSim/models/propulsion/FGPiston.h
@@ -205,8 +205,7 @@ public:
   double getRPM(void) {return RPM;}
 
 protected:
-  double ThrottlePos;
-
+  double ThrottleAngle;
 
 private:
   int crank_counter;
@@ -252,6 +251,7 @@ private:
   //
   double MinManifoldPressure_inHg; // Inches Hg
   double MaxManifoldPressure_inHg; // Inches Hg
+  double MaxManifoldPressure_Percent; // MaxManifoldPressure / 29.92
   double Displacement;             // cubic inches
   double MaxHP;                    // horsepower
   double SparkFailDrop;            // drop of power due to spark failure
@@ -304,7 +304,7 @@ private:
   //
   double rho_air;
   double volumetric_efficiency;
-  double suction_loss;
+  double map_coefficient;
   double m_dot_air;
   double equivalence_ratio;
   double m_dot_fuel;
diff --git a/src/FDM/JSBSim/models/propulsion/FGRocket.cpp b/src/FDM/JSBSim/models/propulsion/FGRocket.cpp
index 66390010c..b4edabfdd 100644
--- a/src/FDM/JSBSim/models/propulsion/FGRocket.cpp
+++ b/src/FDM/JSBSim/models/propulsion/FGRocket.cpp
@@ -106,6 +106,7 @@ FGRocket::~FGRocket(void)
 double FGRocket::Calculate(void)
 {
   double dT = State->Getdt()*Propulsion->GetRate();
+  double thrust;
 
   if (!Flameout && !Starved) ConsumeFuel();
 
@@ -135,7 +136,7 @@ double FGRocket::Calculate(void)
 
     if (Throttle < MinThrottle || Starved) { // Combustion not supported
 
-      PctPower = Thrust = 0.0; // desired thrust
+      PctPower = 0.0; // desired thrust
       Flameout = true;
       VacThrust = 0.0;
 
@@ -149,10 +150,10 @@ double FGRocket::Calculate(void)
 
   } // End thrust calculations
 
-  Thrust = Thruster->Calculate(VacThrust);
-  It += Thrust * dT;
+  thrust = Thruster->Calculate(VacThrust);
+  It += thrust * dT;
 
-  return Thrust;
+  return thrust;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/src/FDM/JSBSim/models/propulsion/FGTank.cpp b/src/FDM/JSBSim/models/propulsion/FGTank.cpp
index caa8ef2f6..af174d46d 100644
--- a/src/FDM/JSBSim/models/propulsion/FGTank.cpp
+++ b/src/FDM/JSBSim/models/propulsion/FGTank.cpp
@@ -148,6 +148,8 @@ FGTank::FGTank(FGFDMExec* exec, Element* el, int tank_number)
   if (Temperature != -9999.0)  InitialTemperature = Temperature = FahrenheitToCelsius(Temperature);
   Area = 40.0 * pow(Capacity/1975, 0.666666667);
 
+  CalculateInertias();
+
   Debug(0);
 }
 
diff --git a/src/FDM/JSBSim/models/propulsion/FGThruster.cpp b/src/FDM/JSBSim/models/propulsion/FGThruster.cpp
index b802d6e8c..9bbe8ed1c 100644
--- a/src/FDM/JSBSim/models/propulsion/FGThruster.cpp
+++ b/src/FDM/JSBSim/models/propulsion/FGThruster.cpp
@@ -86,7 +86,7 @@ FGThruster::FGThruster(FGFDMExec *FDMExec, Element *el, int num ): FGForce(FDMEx
   PropertyManager->Tie( property_name, (FGForce *)this, &FGForce::GetYaw, &FGForce::SetYaw);
 
   if (el->GetName() == "direct") // this is a direct thruster. At this time
-                                      // only a direct thruster can be reversed.
+                                 // only a direct thruster can be reversed.
   {
     snprintf(property_name, 80, "propulsion/engine[%d]/reverser-angle-rad", EngineNum);
     PropertyManager->Tie( property_name, (FGThruster *)this, &FGThruster::GetReverserAngle,
diff --git a/src/FDM/JSBSim/models/propulsion/FGThruster.h b/src/FDM/JSBSim/models/propulsion/FGThruster.h
index c2adf1998..d37a73ad8 100644
--- a/src/FDM/JSBSim/models/propulsion/FGThruster.h
+++ b/src/FDM/JSBSim/models/propulsion/FGThruster.h
@@ -99,7 +99,7 @@ public:
   virtual void SetRPM(double rpm) {};
   virtual double GetPowerRequired(void) {return 0.0;}
   virtual void SetdeltaT(double dt) {deltaT = dt;}
-  double GetThrust(void) {return Thrust;}
+  double GetThrust(void) const {return Thrust;}
   eType GetType(void) {return Type;}
   string GetName(void) {return Name;}
   void SetReverserAngle(double angle) {ReverserAngle = angle;}
diff --git a/src/FDM/JSBSim/models/propulsion/FGTurbine.cpp b/src/FDM/JSBSim/models/propulsion/FGTurbine.cpp
index 022d39235..16bc0fe1e 100644
--- a/src/FDM/JSBSim/models/propulsion/FGTurbine.cpp
+++ b/src/FDM/JSBSim/models/propulsion/FGTurbine.cpp
@@ -107,6 +107,8 @@ void FGTurbine::ResetToIC(void)
 
 double FGTurbine::Calculate(void)
 {
+  double thrust;
+
   TAT = (Auxiliary->GetTotalTemperature() - 491.69) * 0.5555556;
   dt = State->Getdt() * Propulsion->GetRate();
   ThrottlePos = FCS->GetThrottlePos(EngineNumber);
@@ -144,19 +146,19 @@ double FGTurbine::Calculate(void)
   if (Seized) phase = tpSeize;
 
   switch (phase) {
-    case tpOff:    Thrust = Off(); break;
-    case tpRun:    Thrust = Run(); break;
-    case tpSpinUp: Thrust = SpinUp(); break;
-    case tpStart:  Thrust = Start(); break;
-    case tpStall:  Thrust = Stall(); break;
-    case tpSeize:  Thrust = Seize(); break;
-    case tpTrim:   Thrust = Trim(); break;
-    default: Thrust = Off();
+    case tpOff:    thrust = Off(); break;
+    case tpRun:    thrust = Run(); break;
+    case tpSpinUp: thrust = SpinUp(); break;
+    case tpStart:  thrust = Start(); break;
+    case tpStall:  thrust = Stall(); break;
+    case tpSeize:  thrust = Seize(); break;
+    case tpTrim:   thrust = Trim(); break;
+    default: thrust = Off();
   }
 
-  Thrust = Thruster->Calculate(Thrust); // allow thruster to modify thrust (i.e. reversing)
+  thrust = Thruster->Calculate(thrust); // allow thruster to modify thrust (i.e. reversing)
 
-  return Thrust;
+  return thrust;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/src/FDM/JSBSim/models/propulsion/FGTurbine.h b/src/FDM/JSBSim/models/propulsion/FGTurbine.h
index 52193db99..c0524ce1b 100644
--- a/src/FDM/JSBSim/models/propulsion/FGTurbine.h
+++ b/src/FDM/JSBSim/models/propulsion/FGTurbine.h
@@ -169,7 +169,6 @@ public:
   double Calculate(void);
   double CalcFuelNeed(void);
   double GetPowerAvailable(void);
-//  double GetThrust(void) const {return Thrust;}
   /** A lag filter.
       Used to control the rate at which values are allowed to change.
       @param var a pointer to a variable of type double