Sync'ed with JSBSim:

* Fixed the trim on ground algorithm. Now JSBSim aircrafts should no longer be 'dropped' on the runway at start.
* Removed a correction on the propeller induced velocity that was giving erratic results when the aircraft aero velocity is very small.
* Various source comments updates.

src/FDM/JSBSim/initialization/FGTrim.cpp

@@ -57,7 +57,7 @@ using namespace std;
 
 namespace JSBSim {
 
-IDENT(IdSrc,"$Id: FGTrim.cpp,v 1.30 2015/12/29 13:44:39 ehofman Exp $");
+IDENT(IdSrc,"$Id: FGTrim.cpp,v 1.34 2016/06/12 09:09:02 bcoconni Exp $");
 IDENT(IdHdr,ID_TRIM);
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -198,9 +198,8 @@ bool FGTrim::DoTrim(void) {
   double rudder0 = FCS->GetDrCmd();
   double PitchTrim0 = FCS->GetPitchTrimCmd();
 
-  for(int i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++){
+  for(int i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++)
     fdmex->GetGroundReactions()->GetGearUnit(i)->SetReport(false);
-  }
 
   fdmex->SetTrimStatus(true);
   fdmex->SuspendIntegration();
@@ -340,12 +339,12 @@ bool FGTrim::DoTrim(void) {
         cout << endl << "  Trim failed" << endl;
   }
 
+  fdmex->GetPropagate()->InitializeDerivatives();
   fdmex->ResumeIntegration();
   fdmex->SetTrimStatus(false);
 
-  for(int i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++){
+  for(int i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++)
     fdmex->GetGroundReactions()->GetGearUnit(i)->SetReport(true);
-  }
 
   return !trim_failed;
 }
@@ -383,18 +382,22 @@ void FGTrim::trimOnGround(void)
   vector<ContactPoints> contacts;
   FGLocation CGLocation = Propagate->GetLocation();
   FGMatrix33 Tec2b = Propagate->GetTec2b();
-  FGMatrix33 Tl2b = Propagate->GetTl2b();
+  FGMatrix33 Tb2l = Propagate->GetTb2l();
   double hmin = 1E+10;
   int contactRef = -1;
 
   // Build the list of the aircraft contact points and take opportunity of the
   // loop to find which one is closer to (or deeper into) the ground.
-  for (int i = 0; i < GroundReactions->GetNumGearUnits(); i++) {
+  for (int i = 0; i < GroundReactions->GetNumGearUnits(); ++i) {
     ContactPoints c;
     FGLGear* gear = GroundReactions->GetGearUnit(i);
-    c.location = gear->GetLocalGear();
-    FGLocation gearLoc = CGLocation.LocalToLocation(c.location);
-    c.location = Tl2b * c.location;
+
+    // Skip the retracted landing gears
+    if (!gear->GetGearUnitDown())
+      continue;
+
+    c.location = gear->GetBodyLocation();
+    FGLocation gearLoc = CGLocation.LocalToLocation(Tb2l * c.location);
 
     FGColumnVector3 normal, vDummy;
     FGLocation lDummy;
@@ -405,7 +408,7 @@ void FGTrim::trimOnGround(void)
 
     if (height < hmin) {
       hmin = height;
-      contactRef = i;
+      contactRef = contacts.size() - 1;
     }
   }
 
@@ -437,7 +440,7 @@ void FGTrim::trimOnGround(void)
   // Apply the computed rotation to all the contact points
   FGMatrix33 rot = q0.GetTInv();
   vector<ContactPoints>::iterator iter;
-  for (iter = contacts.begin(); iter != contacts.end(); iter++)
+  for (iter = contacts.begin(); iter != contacts.end(); ++iter)
     iter->location = contact0 + rot * (iter->location - contact0);
 
   // Remove the second point to come in contact with the ground from the list.
@@ -461,7 +464,7 @@ void FGTrim::trimOnGround(void)
   FGQuaternion q1(rParam.angleMin, rotationAxis);
 
   // Update the aircraft orientation
-  FGColumnVector3 euler = (q0 * q1 * fgic.GetOrientation()).GetEuler();
+  FGColumnVector3 euler = (fgic.GetOrientation() * q0 * q1).GetEuler();
 
   fgic.SetPhiRadIC(euler(1));
   fgic.SetThetaRadIC(euler(2));
@@ -486,7 +489,7 @@ FGTrim::RotationParameters FGTrim::calcRotation(vector<ContactPoints>& contacts,
 
   rParam.angleMin = 3.0 * M_PI;
 
-  for (iter = contacts.begin(); iter != contacts.end(); iter++) {
+  for (iter = contacts.begin(); iter != contacts.end(); ++iter) {
     // Below the processed contact point is named 'M'
     // Construct an orthonormal basis (u, v, t). The ground normal is obtained
     // from iter->normal.

src/FDM/JSBSim/models/FGAerodynamics.cpp

@@ -50,7 +50,7 @@ using namespace std;
 
 namespace JSBSim {
 
-IDENT(IdSrc,"$Id: FGAerodynamics.cpp,v 1.59 2016/05/23 17:23:36 bcoconni Exp $");
+IDENT(IdSrc,"$Id: FGAerodynamics.cpp,v 1.60 2016/05/30 19:05:58 bcoconni Exp $");
 IDENT(IdHdr,ID_AERODYNAMICS);
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -147,9 +147,13 @@ bool FGAerodynamics::Run(bool Holding)
   unsigned int axis_ctr;
   const double twovel=2*in.Vt;
 
-  // Calculate lift coefficient squared
-  // Make sure that aero/cl-squared is computed with the current qbar
+  // The lift coefficient squared (property aero/cl-squared) is computed before
+  // the aero functions are called to make sure that they use the same value for
+  // qbar.
   if ( in.Qbar > 1.0) {
+    // Skip the computation if qbar is close to zero to avoid huge values for
+    // aero/cl-squared when a non-null lift coincides with a very small aero
+    // velocity (i.e. when qbar is close to zero).
     clsq = (vFw(eLift) + vFwAtCG(eLift))/ (in.Wingarea*in.Qbar);
     clsq *= clsq;
   }

src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.cpp

@@ -153,9 +153,6 @@ double FGStandardAtmosphere::GetPressure(double altitude) const
   double Lmb, Exp, Tmb, deltaH, factor;
   double numRows = StdAtmosTemperatureTable->GetNumRows();
 
-//if (altitude > 328084) // Karman line
-//  altitude = 328084;
-
   // Iterate through the altitudes to find the current Base Altitude
   // in the table. That is, if the current altitude (the argument passed in)
   // is 20000 ft, then the base altitude from the table is 0.0. If the

src/FDM/JSBSim/models/propulsion/FGPropeller.cpp

@@ -45,7 +45,7 @@ using namespace std;
 
 namespace JSBSim {
 
-IDENT(IdSrc,"$Id: FGPropeller.cpp,v 1.57 2016/01/02 17:42:53 bcoconni Exp $");
+IDENT(IdSrc,"$Id: FGPropeller.cpp,v 1.58 2016/06/04 11:06:51 bcoconni Exp $");
 IDENT(IdHdr,ID_PROPELLER);
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -241,20 +241,6 @@ double FGPropeller::Calculate(double EnginePower)
     Vinduced = 0.5 * (-Vel + sqrt(Vel2sum));
   else
     Vinduced = 0.5 * (-Vel - sqrt(-Vel2sum));
-
-  // We need to drop the case where the downstream velocity is opposite in
-  // direction to the aircraft velocity. For example, in such a case, the
-  // direction of the airflow on the tail would be opposite to the airflow on
-  // the wing tips. When such complicated airflows occur, the momentum theory
-  // breaks down and the formulas above are no longer applicable
-  // (see H. Glauert, "The Elements of Airfoil and Airscrew Theory",
-  // 2nd edition, §16.3, pp. 219-221)
-
-  if ((Vel+2.0*Vinduced)*Vel < 0.0) {
-    // The momentum theory is no longer applicable so let's assume the induced
-    // saturates to -0.5*Vel so that the total velocity Vel+2*Vinduced equals 0.
-    Vinduced = -0.5*Vel;
-  }
     
   // P-factor is simulated by a shift of the acting location of the thrust.
   // The shift is a multiple of the angle between the propeller shaft axis