sync with JSB JSBSim CVS

src/FDM/JSBSim/CMakeLists.txt

@@ -77,6 +77,7 @@ set(HEADERS
     models/propulsion/FGRotor.h
     models/propulsion/FGTank.h
     models/propulsion/FGThruster.h
+    models/propulsion/FGTransmission.h
     models/propulsion/FGTurbine.h
     models/propulsion/FGTurboProp.h
     )
@@ -153,6 +154,7 @@ set(SOURCES
     models/propulsion/FGRotor.cpp
     models/propulsion/FGTank.cpp
     models/propulsion/FGThruster.cpp
+    models/propulsion/FGTransmission.cpp
     models/propulsion/FGTurbine.cpp
     models/propulsion/FGTurboProp.cpp
     )

src/FDM/JSBSim/FGFDMExec.cpp

@@ -70,7 +70,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGFDMExec.cpp,v 1.120 2011/11/10 12:06:13 jberndt Exp $";
+static const char *IdSrc = "$Id: FGFDMExec.cpp,v 1.133 2012/04/14 18:10:43 bcoconni Exp $";
 static const char *IdHdr = ID_FDMEXEC;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -98,6 +98,9 @@ FGFDMExec::FGFDMExec(FGPropertyManager* root, unsigned int* fdmctr) : Root(root)
   StandAlone = false;
   firstPass = true;
 
+  IncrementThenHolding = false;  // increment then hold is off by default
+  TimeStepsUntilHold = -1;
+
   sim_time = 0.0;
   dT = 1.0/120.0; // a default timestep size. This is needed for when JSBSim is
                   // run in standalone mode with no initialization file.
@@ -125,7 +128,7 @@ FGFDMExec::FGFDMExec(FGPropertyManager* root, unsigned int* fdmctr) : Root(root)
 
   // Store this FDM's ID
   IdFDM = (*FDMctr); // The main (parent) JSBSim instance is always the "zeroth"
-                                                                      
+
   // Prepare FDMctr for the next child FDM id
   (*FDMctr)++;       // instance. "child" instances are loaded last.
 
@@ -144,9 +147,11 @@ FGFDMExec::FGFDMExec(FGPropertyManager* root, unsigned int* fdmctr) : Root(root)
 
   Constructing = true;
   typedef int (FGFDMExec::*iPMF)(void) const;
+//  typedef unsigned int (FGFDMExec::*uiPMF)(void) const;
 //  instance->Tie("simulation/do_trim_analysis", this, (iPMF)0, &FGFDMExec::DoTrimAnalysis, false);
   instance->Tie("simulation/do_simple_trim", this, (iPMF)0, &FGFDMExec::DoTrim, false);
   instance->Tie("simulation/reset", this, (iPMF)0, &FGFDMExec::ResetToInitialConditions, false);
+  instance->Tie("simulation/randomseed", this, (iPMF)0, &FGFDMExec::SRand, false);
   instance->Tie("simulation/terminate", (int *)&Terminate);
   instance->Tie("simulation/sim-time-sec", this, &FGFDMExec::GetSimTime);
   instance->Tie("simulation/jsbsim-debug", this, &FGFDMExec::GetDebugLevel, &FGFDMExec::SetDebugLevel);
@@ -162,7 +167,7 @@ FGFDMExec::~FGFDMExec()
   try {
     Unbind();
     DeAllocate();
-    
+
     if (IdFDM == 0) { // Meaning this is no child FDM
       if(Root != 0) {
          if(StandAlone)
@@ -408,6 +413,8 @@ void FGFDMExec::LoadInputs(unsigned int idx)
     Propulsion->in.PropFeather      = FCS->GetPropFeather();
     Propulsion->in.H_agl            = Propagate->GetDistanceAGL();
     Propulsion->in.PQR              = Propagate->GetPQR();
+    Propulsion->in.vXYZcg           = MassBalance->GetXYZcg();
+
     break;
   case eAerodynamics:
     Aerodynamics->in.Alpha     = Auxiliary->Getalpha();
@@ -438,9 +445,7 @@ void FGFDMExec::LoadInputs(unsigned int idx)
     GroundReactions->in.TotalDeltaT     = dT * GroundReactions->GetRate();
     GroundReactions->in.WOW             = GroundReactions->GetWOW();
     GroundReactions->in.Location        = Propagate->GetLocation();
-    for (int i=0; i<GroundReactions->GetNumGearUnits(); i++) {
+    GroundReactions->in.vXYZcg          = MassBalance->GetXYZcg();
-      GroundReactions->in.vWhlBodyVec[i] = MassBalance->StructuralToBody(GroundReactions->GetGearUnit(i)->GetLocation());
-    }
     break;
   case eExternalReactions:
     // There are no external inputs to this model.
@@ -479,7 +484,7 @@ void FGFDMExec::LoadInputs(unsigned int idx)
     Accelerations->in.Ti2b     = Propagate->GetTi2b();
     Accelerations->in.Tb2i     = Propagate->GetTb2i();
     Accelerations->in.Tec2b    = Propagate->GetTec2b();
-    Accelerations->in.Tl2b     = Propagate->GetTl2b();
+    Accelerations->in.Tec2i    = Propagate->GetTec2i();
     Accelerations->in.qAttitudeECI = Propagate->GetQuaternionECI();
     Accelerations->in.Moment   = Aircraft->GetMoments();
     Accelerations->in.GroundMoment  = GroundReactions->GetMoments();
@@ -527,9 +532,8 @@ void FGFDMExec::LoadModelConstants(void)
   Aerodynamics->in.Wingspan      = Aircraft->GetWingSpan();
   Auxiliary->in.Wingspan         = Aircraft->GetWingSpan();
   Auxiliary->in.Wingchord        = Aircraft->Getcbar();
-  for (int i=0; i<GroundReactions->GetNumGearUnits(); i++) {
+  GroundReactions->in.vXYZcg     = MassBalance->GetXYZcg();
-    GroundReactions->in.vWhlBodyVec[i] = MassBalance->StructuralToBody(GroundReactions->GetGearUnit(i)->GetLocation());
+  Propulsion->in.vXYZcg          = MassBalance->GetXYZcg();
-  }
 
   LoadPlanetConstants();
 }
@@ -573,10 +577,10 @@ void FGFDMExec::ResetToInitialConditions(int mode)
 {
   if (mode == 1) {
     for (unsigned int i=0; i<Outputs.size(); i++) {
-      Outputs[i]->SetStartNewFile(true); 
+      Outputs[i]->SetStartNewFile(true);
     }
   }
-  
+
   ResetToInitialConditions();
 }
 
@@ -831,21 +835,23 @@ bool FGFDMExec::LoadModel(const string& model, bool addModelToPath)
     while (element) {
       if (debug_lvl > 0) cout << endl << "  Output data set: " << idx << "  ";
       FGOutput* Output = new FGOutput(this);
-      Output->InitModel();
-      Schedule(Output);
       result = Output->Load(element);
       if (!result) {
         cerr << endl << "Aircraft output element has problems in file " << aircraftCfgFileName << endl;
+        delete Output;
         return result;
       } else {
+        Output->InitModel();
+        Schedule(Output);
         Outputs.push_back(Output);
         string outputProp = CreateIndexedPropertyName("simulation/output",idx);
         instance->Tie(outputProp+"/log_rate_hz", Output, (iOPMF)0, &FGOutput::SetRate, false);
-        instance->Tie("simulation/force-output", this, (iOPV)0, &FGFDMExec::ForceOutput, false);
         idx++;
       }
       element = document->FindNextElement("output");
     }
+    if (idx)
+      instance->Tie("simulation/force-output", this, (iOPV)0, &FGFDMExec::ForceOutput, false);
 
     // Lastly, process the child element. This element is OPTIONAL - and NOT YET SUPPORTED.
     element = document->FindElement("child");
@@ -873,7 +879,7 @@ bool FGFDMExec::LoadModel(const string& model, bool addModelToPath)
            << "-------------------------------------------------------------------------------"
            << reset << endl;
     }
-    
+
     if (IsChild) debug_lvl = saved_debug_lvl;
 
   } else {
@@ -1072,7 +1078,7 @@ bool FGFDMExec::ReadChild(Element* el)
     cerr << endl << highint << fgred << "  No location was found for this child object!" << reset << endl;
     exit(-1);
   }
-  
+
   Element* orientation = el->FindElement("orient");
   if (orientation) {
     child->Orient = orientation->FindElementTripletConvertTo("RAD");
@@ -1121,6 +1127,30 @@ void FGFDMExec::EnableOutput(void)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+void FGFDMExec::CheckIncrementalHold(void)
+{
+  // Only check if increment then hold is on
+  if( IncrementThenHolding ) {
+
+    if (TimeStepsUntilHold == 0) {
+
+      // Should hold simulation if TimeStepsUntilHold has reached zero
+      holding = true;
+
+      // Still need to decrement TimeStepsUntilHold as value of -1
+      // indicates that incremental then hold is turned off
+      IncrementThenHolding = false;
+      TimeStepsUntilHold--;
+
+    } else if ( TimeStepsUntilHold > 0 ) {
+      // Keep decrementing until 0 is reached
+      TimeStepsUntilHold--;
+    }
+  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 void FGFDMExec::ForceOutput(int idx)
 {
   if (idx >= (int)0 && idx < (int)Outputs.size()) Outputs[idx]->Print();
@@ -1134,11 +1164,11 @@ bool FGFDMExec::SetOutputDirectives(const string& fname)
 
   FGOutput* Output = new FGOutput(this);
   Output->SetDirectivesFile(RootDir + fname);
-  Output->InitModel();
-  Schedule(Output);
   result = Output->Load(0);
 
   if (result) {
+    Output->InitModel();
+    Schedule(Output);
     Output->Run(holding);
     Outputs.push_back(Output);
     typedef double (FGOutput::*iOPMF)(void) const;
@@ -1170,6 +1200,13 @@ void FGFDMExec::DoTrim(int mode)
   sim_time = saved_time;
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void FGFDMExec::SRand(int sr)
+{
+  srand(sr);
+}
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 //    The bitmasked value choices are as follows:
 //    unset: In this case (the default) JSBSim would only print
@@ -1195,7 +1232,7 @@ void FGFDMExec::Debug(int from)
 
   if (debug_lvl & 1 && IdFDM == 0) { // Standard console startup message output
     if (from == 0) { // Constructor
-      cout << "\n\n     " 
+      cout << "\n\n     "
            << "JSBSim Flight Dynamics Model v" << JSBSim_version << endl;
       cout << "            [JSBSim-ML v" << needed_cfg_version << "]\n\n";
       cout << "JSBSim startup beginning ...\n\n";

src/FDM/JSBSim/FGFDMExec.h

@@ -55,7 +55,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_FDMEXEC "$Id: FGFDMExec.h,v 1.74 2011/11/09 21:58:26 bcoconni Exp $"
+#define ID_FDMEXEC "$Id: FGFDMExec.h,v 1.76 2012/01/21 16:46:08 jberndt Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -180,7 +180,7 @@ CLASS DOCUMENTATION
                                 property actually maps toa function call of DoTrim().
 
     @author Jon S. Berndt
-    @version $Revision: 1.74 $
+    @version $Revision: 1.76 $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -464,6 +464,10 @@ public:
   void EnableOutput(void);
   /// Pauses execution by preventing time from incrementing.
   void Hold(void) {holding = true;}
+  /// Turn on hold after increment
+  void EnableIncrementThenHold(int Timesteps) {TimeStepsUntilHold = Timesteps; IncrementThenHolding = true;}
+  /// Checks if required to hold afer increment
+  void CheckIncrementalHold(void);
   /// Resumes execution from a "Hold".
   void Resume(void) {holding = false;}
   /// Returns true if the simulation is Holding (i.e. simulation time is not moving).
@@ -563,6 +567,8 @@ private:
   double saved_dT;
   double sim_time;
   bool holding;
+  bool IncrementThenHolding;
+  int TimeStepsUntilHold;
   bool Constructing;
   bool modelLoaded;
   bool IsChild;
@@ -615,6 +621,7 @@ private:
   bool ReadChild(Element*);
   bool ReadPrologue(Element*);
   void ResetToInitialConditions(int mode);
+  void SRand(int sr);
   void LoadInputs(unsigned int idx);
   void LoadPlanetConstants(void);
   void LoadModelConstants(void);

src/FDM/JSBSim/FGJSBBase.cpp

@@ -44,7 +44,7 @@ INCLUDES
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGJSBBase.cpp,v 1.32 2011/10/22 14:38:30 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGJSBBase.cpp,v 1.35 2012/03/25 11:05:36 bcoconni Exp $";
 static const char *IdHdr = ID_JSBBASE;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -159,7 +159,7 @@ void FGJSBBase::PutMessage(const string& text, int iVal)
   msg.messageId = messageId++;
   msg.subsystem = "FDM";
   msg.type = Message::eInteger;
-  msg.bVal = (iVal != 0);
+  msg.iVal = iVal;
   Messages.push(msg);
 }
 
@@ -172,7 +172,7 @@ void FGJSBBase::PutMessage(const string& text, double dVal)
   msg.messageId = messageId++;
   msg.subsystem = "FDM";
   msg.type = Message::eDouble;
-  msg.bVal = (dVal != 0.0);
+  msg.dVal = dVal;
   Messages.push(msg);
 }
 

src/FDM/JSBSim/JSBSim.cxx

@@ -378,7 +378,7 @@ void FGJSBsim::init()
 
     if (fgGetBool("/environment/params/control-fdm-atmosphere")) {
       Atmosphere->SetTemperature(temperature->getDoubleValue(), get_Altitude(), FGAtmosphere::eCelsius);
-      Atmosphere->SetPressureSL(pressureSL->getDoubleValue(), FGAtmosphere::eInchesHg);
+      Atmosphere->SetPressureSL(FGAtmosphere::eInchesHg, pressureSL->getDoubleValue());
       // initialize to no turbulence, these values get set in the update loop
       Winds->SetTurbType(FGWinds::ttNone);
       Winds->SetTurbGain(0.0);
@@ -666,7 +666,7 @@ bool FGJSBsim::copy_to_JSBsim()
     }
 
     Atmosphere->SetTemperature(temperature->getDoubleValue(), get_Altitude(), FGAtmosphere::eCelsius);
-    Atmosphere->SetPressureSL(pressureSL->getDoubleValue(), FGAtmosphere::eInchesHg);
+    Atmosphere->SetPressureSL(FGAtmosphere::eInchesHg, pressureSL->getDoubleValue());
 
     Winds->SetTurbType((FGWinds::tType)TURBULENCE_TYPE_NAMES[turbulence_model->getStringValue()]);
     switch( Winds->GetTurbType() ) {
@@ -1222,7 +1222,7 @@ void FGJSBsim::init_gear(void )
       node->setDoubleValue("rollspeed-ms", gear->GetWheelRollVel()*0.3043);
       node->setBoolValue("has-brake", gear->GetBrakeGroup() > 0);
       node->setDoubleValue("position-norm", gear->GetGearUnitPos());
-      node->setDoubleValue("tire-pressure-norm", gear->GetTirePressure());
+//    node->setDoubleValue("tire-pressure-norm", gear->GetTirePressure());
       node->setDoubleValue("compression-norm", gear->GetCompLen());
       node->setDoubleValue("compression-ft", gear->GetCompLen());
       if ( gear->GetSteerable() )
@@ -1240,7 +1240,7 @@ void FGJSBsim::update_gear(void)
       node->getChild("wow", 0, true)->setBoolValue( gear->GetWOW());
       node->getChild("rollspeed-ms", 0, true)->setDoubleValue(gear->GetWheelRollVel()*0.3043);
       node->getChild("position-norm", 0, true)->setDoubleValue(gear->GetGearUnitPos());
-      gear->SetTirePressure(node->getDoubleValue("tire-pressure-norm"));
+//    gear->SetTirePressure(node->getDoubleValue("tire-pressure-norm"));
       node->setDoubleValue("compression-norm", gear->GetCompLen());
       node->setDoubleValue("compression-ft", gear->GetCompLen());
       if ( gear->GetSteerable() )

src/FDM/JSBSim/initialization/FGInitialCondition.cpp

@@ -63,7 +63,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGInitialCondition.cpp,v 1.78 2011/11/09 21:57:51 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGInitialCondition.cpp,v 1.81 2012/04/08 15:22:56 jberndt Exp $";
 static const char *IdHdr = ID_INITIALCONDITION;
 
 //******************************************************************************
@@ -165,10 +165,10 @@ void FGInitialCondition::WriteStateFile(int num)
     filename = "initfile.xml";
   else
     filename.append("/initfile.xml");
-  
+
   ofstream outfile(filename.c_str());
   FGPropagate* Propagate = fdmex->GetPropagate();
-  
+
   if (outfile.is_open()) {
     outfile << "<?xml version=\"1.0\"?>" << endl;
     outfile << "<initialize name=\"reset00\">" << endl;
@@ -1026,13 +1026,14 @@ bool FGInitialCondition::Load_v1(void)
 
   // Refer to Stevens and Lewis, 1.5-14a, pg. 49.
   // This is the rotation rate of the "Local" frame, expressed in the local frame.
+  const FGMatrix33& Tl2b = orientation.GetT();
   double radInv = 1.0 / position.GetRadius();
   FGColumnVector3 vOmegaLocal = FGColumnVector3(
    radInv*vUVW_NED(eEast),
   -radInv*vUVW_NED(eNorth),
   -radInv*vUVW_NED(eEast)*position.GetTanLatitude() );
 
-  vPQR_body = vOmegaLocal;
+  vPQR_body = Tl2b * vOmegaLocal;
 
   return result;
 }
@@ -1182,7 +1183,7 @@ bool FGInitialCondition::Load_v2(void)
   // - Local
   // - Body
   // The vehicle will be defaulted to (0,0,0) in the Body frame if nothing is provided.
-  
+
   Element* velocity_el = document->FindElement("velocity");
   FGColumnVector3 vInitVelocity = FGColumnVector3(0.0, 0.0, 0.0);
   FGMatrix33 mTec2l = position.GetTec2l();
@@ -1230,7 +1231,7 @@ bool FGInitialCondition::Load_v2(void)
   // - ECI (Earth Centered Inertial)
   // - ECEF (Earth Centered, Earth Fixed)
   // - Body
-  
+
   Element* attrate_el = document->FindElement("attitude_rate");
   const FGMatrix33& Tl2b = orientation.GetT();
 
@@ -1253,19 +1254,21 @@ bool FGInitialCondition::Load_v2(void)
     } else if (frame == "ecef") {
       vPQR_body = Tl2b * position.GetTec2l() * vAttRate;
     } else if (frame == "local") {
-      vPQR_body = vAttRate + vOmegaLocal;
+      vPQR_body = Tl2b * (vAttRate + vOmegaLocal);
+    } else if (frame == "body") {
+      vPQR_body = vAttRate;
     } else if (!frame.empty()) { // misspelling of frame
-      
+
       cerr << endl << fgred << "  Attitude rate frame type: \"" << frame
            << "\" is not supported!" << reset << endl << endl;
       result = false;
 
     } else if (frame.empty()) {
-      vPQR_body = vOmegaLocal;
+      vPQR_body = Tl2b * vOmegaLocal;
     }
 
   } else { // Body frame attitude rate assumed 0 relative to local.
-      vPQR_body = vOmegaLocal;
+      vPQR_body = Tl2b * vOmegaLocal;
   }
 
   return result;

src/FDM/JSBSim/initialization/FGInitialCondition.h

@@ -49,12 +49,13 @@ INCLUDES
 
 #include "input_output/FGXMLFileRead.h"
 #include "math/FGLocation.h"
+#include "math/FGQuaternion.h"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_INITIALCONDITION "$Id: FGInitialCondition.h,v 1.32 2011/11/06 18:14:51 bcoconni Exp $"
+#define ID_INITIALCONDITION "$Id: FGInitialCondition.h,v 1.33 2011/12/10 15:49:21 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -215,7 +216,7 @@ CLASS DOCUMENTATION
    @property ic/r-rad_sec (read/write) Yaw rate initial condition in radians/second
 
    @author Tony Peden
-   @version "$Id: FGInitialCondition.h,v 1.32 2011/11/06 18:14:51 bcoconni Exp $"
+   @version "$Id: FGInitialCondition.h,v 1.33 2011/12/10 15:49:21 bcoconni Exp $"
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/FDM/JSBSim/math/FGColumnVector3.cpp

@@ -47,7 +47,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGColumnVector3.cpp,v 1.14 2010/12/07 12:57:14 jberndt Exp $";
+static const char *IdSrc = "$Id: FGColumnVector3.cpp,v 1.15 2012/02/07 00:27:51 jentron Exp $";
 static const char *IdHdr = ID_COLUMNVECTOR3;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -65,9 +65,9 @@ FGColumnVector3::FGColumnVector3(void)
 string FGColumnVector3::Dump(const string& delimiter) const
 {
   ostringstream buffer;
-  buffer << std::setw(18) << std::setprecision(16) << data[0] << delimiter;
+  buffer << std::setprecision(16) << data[0] << delimiter;
-  buffer << std::setw(18) << std::setprecision(16) << data[1] << delimiter;
+  buffer << std::setprecision(16) << data[1] << delimiter;
-  buffer << std::setw(18) << std::setprecision(16) << data[2];
+  buffer << std::setprecision(16) << data[2];
   return buffer.str();
 }
 

src/FDM/JSBSim/math/FGColumnVector3.h

@@ -46,7 +46,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_COLUMNVECTOR3 "$Id: FGColumnVector3.h,v 1.16 2010/12/07 12:57:14 jberndt Exp $"
+#define ID_COLUMNVECTOR3 "$Id: FGColumnVector3.h,v 1.17 2011/12/10 15:49:21 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -60,7 +60,7 @@ CLASS DOCUMENTATION
 
 /** This class implements a 3 element column vector.
     @author Jon S. Berndt, Tony Peden, et. al.
-    @version $Id: FGColumnVector3.h,v 1.16 2010/12/07 12:57:14 jberndt Exp $
+    @version $Id: FGColumnVector3.h,v 1.17 2011/12/10 15:49:21 bcoconni Exp $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -243,7 +243,7 @@ public:
   /** Dot product of two vectors
       Compute and return the euclidean dot (or scalar) product of two vectors
       v1 and v2 */
-  friend inline double DotProduct(const FGColumnVector3& v1, const FGColumnVector3& v2) {
+  friend double DotProduct(const FGColumnVector3& v1, const FGColumnVector3& v2) {
     return v1.data[0]*v2.data[0] + v1.data[1]*v2.data[1] + v1.data[2]*v2.data[2];
   }
 

src/FDM/JSBSim/math/FGFunction.cpp

@@ -43,7 +43,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGFunction.cpp,v 1.42 2011/09/07 02:36:04 jberndt Exp $";
+static const char *IdSrc = "$Id: FGFunction.cpp,v 1.43 2012/02/05 11:15:54 bcoconni Exp $";
 static const char *IdHdr = ID_FUNCTION;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -476,7 +476,7 @@ double FGFunction::GetValue(void) const
   case eIfThen:
     {
       i = Parameters.size();
-      if (i != 3) {
+      if (i == 3) {
         if (GetBinary(temp) == 1) {
           temp = Parameters[1]->GetValue();
         } else {

src/FDM/JSBSim/math/FGLocation.cpp

@@ -48,7 +48,7 @@ INCLUDES
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGLocation.cpp,v 1.26 2011/11/06 18:14:51 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGLocation.cpp,v 1.29 2012/04/14 12:14:37 bcoconni Exp $";
 static const char *IdHdr = ID_LOCATION;
 using std::cerr;
 using std::endl;
@@ -61,9 +61,8 @@ CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 FGLocation::FGLocation(void)
+  : mECLoc(1.0, 0.0, 0.0), mCacheValid(false)
 {
-  mCacheValid = false;
-
   a = b = a2 = b2 = 0.0;
   e = e2 = f = 1.0;
   eps2 = -1.0;
@@ -78,15 +77,13 @@ FGLocation::FGLocation(void)
   mTec2i.InitMatrix();
   mTi2l.InitMatrix();
   mTl2i.InitMatrix();
-  mECLoc.InitMatrix();
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 FGLocation::FGLocation(double lon, double lat, double radius)
+  : mCacheValid(false)
 {
-  mCacheValid = false;
-
   a = b = a2 = b2 = 0.0;
   e = e2 = f = 1.0;
   eps2 = -1.0;
@@ -113,7 +110,8 @@ FGLocation::FGLocation(double lon, double lat, double radius)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-FGLocation::FGLocation(const FGColumnVector3& lv) : mECLoc(lv), mCacheValid(false)
+FGLocation::FGLocation(const FGColumnVector3& lv)
+  : mECLoc(lv), mCacheValid(false)
 {
   a = b = a2 = b2 = 0.0;
   e = e2 = f = 1.0;
@@ -136,7 +134,6 @@ FGLocation::FGLocation(const FGColumnVector3& lv) : mECLoc(lv), mCacheValid(fals
 FGLocation::FGLocation(const FGLocation& l)
   : mECLoc(l.mECLoc), mCacheValid(l.mCacheValid)
 {
-
   a = l.a;
   b = l.b;
   a2 = l.a2;
@@ -289,7 +286,7 @@ void FGLocation::SetPosition(double lon, double lat, double radius)
 void FGLocation::SetPositionGeodetic(double lon, double lat, double height)
 {
   mCacheValid = false;
-  
+
   mGeodLat = lat;
   mLon = lon;
   GeodeticAltitude = height;
@@ -319,30 +316,6 @@ void FGLocation::SetEllipse(double semimajor, double semiminor)
   f = 1.0 - b/a;
 }
 
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-// Compute the ECEF to ECI transformation matrix using Stevens and Lewis "Aircraft
-// Control and Simulation", second edition, eqn. 1.4-12, pg. 39. In Stevens and Lewis
-// notation, this is C_i/e, a transformation from ECEF to ECI.
-
-const FGMatrix33& FGLocation::GetTec2i(void)
-{
-  ComputeDerived();
-  return mTec2i;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-// This is given in Stevens and Lewis "Aircraft
-// Control and Simulation", second edition, eqn. 1.4-12, pg. 39
-// The notation in Stevens and Lewis is: C_e/i. This represents a transformation
-// from ECI to ECEF - and the orientation of the ECEF frame relative to the ECI
-// frame.
-
-const FGMatrix33& FGLocation::GetTi2ec(void)
-{
-  ComputeDerived();
-  return mTi2ec;
-}
-
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 void FGLocation::ComputeDerivedUnconditional(void) const
@@ -388,7 +361,7 @@ void FGLocation::ComputeDerivedUnconditional(void) const
 
   // Compute the transform matrices from and to the earth centered frame.
   // See Stevens and Lewis, "Aircraft Control and Simulation", Second Edition,
-  // Eqn. 1.4-13, page 40. In Stevens and Lewis notation, this is C_n/e - the 
+  // Eqn. 1.4-13, page 40. In Stevens and Lewis notation, this is C_n/e - the
   // orientation of the navigation (local) frame relative to the ECEF frame,
   // and a transformation from ECEF to nav (local) frame.
 
@@ -396,7 +369,7 @@ void FGLocation::ComputeDerivedUnconditional(void) const
                            -sinLon   ,     cosLon    ,    0.0 ,
                        -cosLon*cosLat, -sinLon*cosLat, -sinLat  );
 
-  // In Stevens and Lewis notation, this is C_e/n - the 
+  // In Stevens and Lewis notation, this is C_e/n - the
   // orientation of the ECEF frame relative to the nav (local) frame,
   // and a transformation from nav (local) to ECEF frame.
 
@@ -422,7 +395,7 @@ void FGLocation::ComputeDerivedUnconditional(void) const
 
   if (a != 0.0 && b != 0.0) {
     double c, p, q, s, t, u, v, w, z, p2, u2, r0;
-    double Ne, P, Q0, Q, signz0, sqrt_q, z_term; 
+    double Ne, P, Q0, Q, signz0, sqrt_q, z_term;
     p  = fabs(mECLoc(eZ))/eps2;
     s  = r02/(e2*eps2);
     p2 = p*p;

src/FDM/JSBSim/math/FGLocation.h

@@ -52,7 +52,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_LOCATION "$Id: FGLocation.h,v 1.29 2011/11/06 18:14:51 bcoconni Exp $"
+#define ID_LOCATION "$Id: FGLocation.h,v 1.31 2012/02/05 14:56:17 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -65,7 +65,7 @@ CLASS DOCUMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 /** FGLocation holds an arbitrary location in the Earth centered Earth fixed
-    reference frame (ECEF). This coordinate frame has its center in the middle
+    reference frame (ECEF). The coordinate frame ECEF has its center in the middle
     of the earth. The X-axis points from the center of the Earth towards a
     location with zero latitude and longitude on the Earth surface. The Y-axis
     points from the center of the Earth towards a location with zero latitude
@@ -95,10 +95,15 @@ CLASS DOCUMENTATION
     conversion functions for conversion of position vectors given in the one
     frame to positions in the other frame.
 
+    To keep the transformation matrices between the ECI and ECEF frames up to
+    date, the Earth angular position must be updated by calling
+    SetEarthPositionAngle() or IncrementEarthPositionAngle(). This must be done
+    prior to any conversion from and to the ECI frame.
+
     The Earth centered reference frame is NOT an inertial frame since it rotates
     with the Earth.
 
-    The coordinates in the Earth centered frame are the master values. All other
+    The cartesian coordinates (X,Y,Z) in the Earth centered frame are the master values. All other
     values are computed from these master values and are cached as long as the
     location is changed by access through a non-const member function. Values
     are cached to improve performance. It is best practice to work with a
@@ -112,14 +117,14 @@ CLASS DOCUMENTATION
 
     Given,
 
-    -that we model a vehicle near the Earth
+    - that we model a vehicle near the Earth
-    -that the Earth surface radius is about 2*10^7, ft
+    - that the Earth surface radius is about 2*10^7, ft
-    -that we use double values for the representation of the location
+    - that we use double values for the representation of the location
-    
+
     we have an accuracy of about
-    
+
     1e-16*2e7ft/1 = 2e-9 ft
-    
+
     left. This should be sufficient for our needs. Note that this is the same
     relative accuracy we would have when we compute directly with
     lon/lat/radius. For the radius value this is clear. For the lon/lat pair
@@ -146,7 +151,7 @@ CLASS DOCUMENTATION
     @see W. C. Durham "Aircraft Dynamics & Control", section 2.2
 
     @author Mathias Froehlich
-    @version $Id: FGLocation.h,v 1.29 2011/11/06 18:14:51 bcoconni Exp $
+    @version $Id: FGLocation.h,v 1.31 2012/02/05 14:56:17 bcoconni Exp $
   */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -166,7 +171,10 @@ public:
       @param radius distance from center of earth to vehicle in feet*/
   FGLocation(double lon, double lat, double radius);
 
-  /** Column constructor. */
+  /** Constructor to initialize the location with the cartesian coordinates
+      (X,Y,Z) contained in the input FGColumnVector3. Distances are in feet,
+      the position is expressed in the ECEF frame.
+      @param lv vector that contain the cartesian coordinates*/
   FGLocation(const FGColumnVector3& lv);
 
   /** Copy constructor. */
@@ -301,6 +309,11 @@ public:
       return -mTec2l(3,3)/cLat;
   }
 
+  /** Return the Earth Position Angle.
+      This is the relative orientation of the ECEF frame with respect to the
+      Inertial frame.
+      @return the Earth fixed frame (ECEF) rotation offset about the axis with
+              respect to the Inertial (ECI) frame in radians. */
   double GetEPA() const {return epa;}
 
   /** Get the distance from the center of the earth.
@@ -310,35 +323,46 @@ public:
   //double GetRadius() const { return mECLoc.Magnitude(); } // may not work with FlightGear
   double GetRadius() const { ComputeDerived(); return mRadius; }
 
-  /// @name Functions that need the ground callback to be set
+  /** @name Functions that rely on the ground callback
+      The following functions allow to set and get the vehicle position above
+      the sea or the ground. The sea and the ground levels are obtained by
+      interrogating an FGGroundCallback instance. A ground callback must
+      therefore be set with SetGroundCallback() before calling any of these
+      functions. */
   ///@{
   /** Set the altitude above sea level.
-      @param altitudeASL altitude above Sea Level in feet. */
+      @param altitudeASL altitude above Sea Level in feet.
+      @see SetGroundCallback */
   void SetAltitudeASL(double altitudeASL)
   { SetRadius(GroundCallback->GetSeaLevelRadius(*this) + altitudeASL); }
 
   /** Set the altitude above ground level.
-      @param altitudeAGL altitude above Ground Level in feet. */
+      @param altitudeAGL altitude above Ground Level in feet.
+      @see SetGroundCallback */
   void SetAltitudeAGL(double altitudeAGL, double time)
   { SetRadius(GroundCallback->GetTerrainGeoCentRadius(time, *this) + altitudeAGL); }
 
   /** Get the local sea level radius
-      @return the sea level radius at the location in feet. */
+      @return the sea level radius at the location in feet.
+      @see SetGroundCallback */
   double GetSeaLevelRadius(void) const
   { ComputeDerived(); return GroundCallback->GetSeaLevelRadius(*this); }
 
   /** Get the local terrain radius
-      @return the terrain level radius at the location in feet. */
+      @return the terrain level radius at the location in feet.
+      @see SetGroundCallback */
   double GetTerrainRadius(double time) const
   { ComputeDerived(); return GroundCallback->GetTerrainGeoCentRadius(time, *this); }
 
   /** Get the altitude above sea level.
-      @return the altitude ASL in feet. */
+      @return the altitude ASL in feet.
+      @see SetGroundCallback */
   double GetAltitudeASL() const
   { ComputeDerived(); return GroundCallback->GetAltitude(*this); }
 
   /** Get the altitude above ground level.
-      @return the altitude AGL in feet. */
+      @return the altitude AGL in feet.
+      @see SetGroundCallback */
   double GetAltitudeAGL(double time) const {
     FGLocation c;
     FGColumnVector3 n,v,w;
@@ -351,55 +375,79 @@ public:
       @param normal  Terrain normal vector in contact point    (ECEF frame)
       @param v       Terrain linear velocity in contact point  (ECEF frame)
       @param w       Terrain angular velocity in contact point (ECEF frame)
-      @return Location altitude above contact point (AGL) in feet. */
+      @return Location altitude above contact point (AGL) in feet.
+      @see SetGroundCallback */
   double GetContactPoint(double time,
                          FGLocation& contact, FGColumnVector3& normal,
                          FGColumnVector3& v, FGColumnVector3& w) const
   { ComputeDerived(); return GroundCallback->GetAGLevel(time, *this, contact, normal, v, w); }
+  ///@}
 
-  /** Sets the ground callback pointer. For optimal memory management, a shared
+  /** Sets the ground callback pointer. The FGGroundCallback instance will be
-      pointer is used internally that maintains a reference counter. The calling
+      interrogated by FGLocation each time some terrain informations are needed.
-      application must therefore use FGGroundCallback_ptr 'smart pointers' to
+      This will mainly occur when altitudes above the sea level or above the
-      manage their copy of the ground callback.
+      ground level are needed. A 'smart pointer' is used internally to prevent
+      the FGGroundCallback instance against accidental deletion. This can only
+      work if the calling application also make use of FGGroundCallback_ptr
+      'smart pointers' to manage their copy of the ground callback.
       @param gc A pointer to a ground callback object
       @see FGGroundCallback
    */
   static void SetGroundCallback(FGGroundCallback* gc) { GroundCallback = gc; }
 
-  /** Get a pointer to the ground callback currently used. It is recommanded
+  /** Get a pointer to the ground callback currently used. Since the
-      to store the returned pointer in a 'smart pointer' FGGroundCallback_ptr.
+      FGGroundcallback instance might have been created outside JSBSim, it is
+      recommanded to store the returned pointer in a 'smart pointer'
+      FGGroundCallback_ptr. This pointer maintains a reference counter and
+      protects the returned pointer against an accidental deletion of the object
+      it is pointing to.
       @return A pointer to the current ground callback object.
       @see FGGroundCallback
    */
   static FGGroundCallback* GetGroundCallback(void) { return GroundCallback; }
-  ///@}
 
   /** Transform matrix from local horizontal to earth centered frame.
-      Returns a const reference to the rotation matrix of the transform from
+      @return a const reference to the rotation matrix of the transform from
       the local horizontal frame to the earth centered frame. */
   const FGMatrix33& GetTl2ec(void) const { ComputeDerived(); return mTl2ec; }
 
   /** Transform matrix from the earth centered to local horizontal frame.
-      Returns a const reference to the rotation matrix of the transform from
+      @return a const reference to the rotation matrix of the transform from
       the earth centered frame to the local horizontal frame. */
   const FGMatrix33& GetTec2l(void) const { ComputeDerived(); return mTec2l; }
 
   /** Transform matrix from inertial to earth centered frame.
-      Returns a const reference to the rotation matrix of the transform from
+      @return a const reference to the rotation matrix of the transform from
-      the inertial frame to the earth centered frame (ECI to ECEF). */
+      the inertial frame to the earth centered frame (ECI to ECEF).
-  const FGMatrix33& GetTi2ec(void);
+      @see SetEarthPositionAngle
+      @see IncrementEarthPositionAngle */
+  const FGMatrix33& GetTi2ec(void) const { ComputeDerived(); return mTi2ec; }
 
   /** Transform matrix from the earth centered to inertial frame.
-      Returns a const reference to the rotation matrix of the transform from
+      @return a const reference to the rotation matrix of the transform from
-      the earth centered frame to the inertial frame (ECEF to ECI). */
+      the earth centered frame to the inertial frame (ECEF to ECI).
-  const FGMatrix33& GetTec2i(void);
+      @see SetEarthPositionAngle
+      @see IncrementEarthPositionAngle */
+  const FGMatrix33& GetTec2i(void) const { ComputeDerived(); return mTec2i; }
 
+  /** Transform matrix from the inertial to local horizontal frame.
+      @return a const reference to the rotation matrix of the transform from
+      the inertial frame to the local horizontal frame.
+      @see SetEarthPositionAngle
+      @see IncrementEarthPositionAngle */
   const FGMatrix33& GetTi2l(void) const {ComputeDerived(); return mTi2l;}
 
+  /** Transform matrix from local horizontal to inertial frame.
+      @return a const reference to the rotation matrix of the transform from
+      the local horizontal frame to the inertial frame.
+      @see SetEarthPositionAngle
+      @see IncrementEarthPositionAngle */
   const FGMatrix33& GetTl2i(void) const {ComputeDerived(); return mTl2i;}
 
   /** Conversion from Local frame coordinates to a location in the
       earth centered and fixed frame.
+      This function calculates the FGLocation of an object which position
+      relative to the vehicle is given as in input.
       @param lvec Vector in the local horizontal coordinate frame
       @return The location in the earth centered and fixed frame */
   FGLocation LocalToLocation(const FGColumnVector3& lvec) const {
@@ -408,6 +456,8 @@ public:
 
   /** Conversion from a location in the earth centered and fixed frame
       to local horizontal frame coordinates.
+      This function calculates the relative position between the vehicle and
+      the input vector and returns the result expressed in the local frame.
       @param ecvec Vector in the earth centered and fixed frame
       @return The vector in the local horizontal coordinate frame */
   FGColumnVector3 LocationToLocal(const FGColumnVector3& ecvec) const {
@@ -456,7 +506,7 @@ public:
       The location can be set by an Earth-centered, Earth-fixed (ECEF) frame
       position vector. The cache is marked as invalid, so any future requests
       for selected important data will cause the parameters to be calculated.
-      @param v the ECEF column vector in feet. 
+      @param v the ECEF column vector in feet.
       @return a reference to the FGLocation object. */
   const FGLocation& operator=(const FGColumnVector3& v)
   {
@@ -464,12 +514,12 @@ public:
     mECLoc(eY) = v(eY);
     mECLoc(eZ) = v(eZ);
     mCacheValid = false;
-    ComputeDerived();
+    //ComputeDerived();
     return *this;
   }
 
   /** Sets this location via the supplied location object.
-      @param v A location object reference. 
+      @param v A location object reference.
       @return a reference to the FGLocation object. */
   const FGLocation& operator=(const FGLocation& l);
 
@@ -484,39 +534,61 @@ public:
   bool operator!=(const FGLocation& l) const { return ! operator==(l); }
 
   /** This operator adds the ECEF position vectors.
-      The supplied vector (right side) is added to the ECEF position vector
+      The cartesian coordinates of the supplied vector (right side) are added to
-      on the left side of the equality, and a pointer to this object is
+      the ECEF position vector on the left side of the equality, and a reference
-      returned. */
+      to this object is returned. */
   const FGLocation& operator+=(const FGLocation &l) {
     mCacheValid = false;
     mECLoc += l.mECLoc;
     return *this;
   }
 
+  /** This operator substracts the ECEF position vectors.
+      The cartesian coordinates of the supplied vector (right side) are
+      substracted from the ECEF position vector on the left side of the
+      equality, and a reference to this object is returned. */
   const FGLocation& operator-=(const FGLocation &l) {
     mCacheValid = false;
     mECLoc -= l.mECLoc;
     return *this;
   }
 
+  /** This operator scales the ECEF position vector.
+      The cartesian coordinates of the ECEF position vector on the left side of
+      the equality are scaled by the supplied value (right side), and a
+      reference to this object is returned. */
   const FGLocation& operator*=(double scalar) {
     mCacheValid = false;
     mECLoc *= scalar;
     return *this;
   }
 
+  /** This operator scales the ECEF position vector.
+      The cartesian coordinates of the ECEF position vector on the left side of
+      the equality are scaled by the inverse of the supplied value (right side),
+      and a reference to this object is returned. */
   const FGLocation& operator/=(double scalar) {
     return operator*=(1.0/scalar);
   }
 
+  /** This operator adds two ECEF position vectors.
+      A new object is returned that defines a position which is the sum of the
+      cartesian coordinates of the two positions provided. */
   FGLocation operator+(const FGLocation& l) const {
     return FGLocation(mECLoc + l.mECLoc);
   }
 
+  /** This operator substracts two ECEF position vectors.
+      A new object is returned that defines a position which is the difference
+      of the cartesian coordinates of the two positions provided. */
   FGLocation operator-(const FGLocation& l) const {
     return FGLocation(mECLoc - l.mECLoc);
   }
 
+  /** This operator scales an ECEF position vector.
+      A new object is returned that defines a position made of the cartesian
+      coordinates of the provided ECEF position scaled by the supplied scalar
+      value. */
   FGLocation operator*(double scalar) const {
     return FGLocation(scalar*mECLoc);
   }
@@ -558,7 +630,7 @@ private:
   mutable double mRadius;
   mutable double mGeodLat;
   mutable double GeodeticAltitude;
-  
+
   double initial_longitude;
 
   /** The cached rotation matrices from and to the associated frames. */

src/FDM/JSBSim/math/FGMatrix33.cpp

@@ -39,6 +39,7 @@ INCLUDES
 
 #include "FGMatrix33.h"
 #include "FGColumnVector3.h"
+#include "FGQuaternion.h"
 #include <sstream>
 #include <iomanip>
 
@@ -48,7 +49,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGMatrix33.cpp,v 1.11 2010/12/07 12:57:14 jberndt Exp $";
+static const char *IdSrc = "$Id: FGMatrix33.cpp,v 1.12 2011/12/10 15:49:21 bcoconni Exp $";
 static const char *IdHdr = ID_MATRIX33;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -119,7 +120,7 @@ FGQuaternion FGMatrix33::GetQuaternion(void)
 
   // Find largest of the above
   idx = 0;
-  for (int i=1; i<4; i++) if (tempQ[i] > tempQ[idx]) idx = i; 
+  for (int i=1; i<4; i++) if (tempQ[i] > tempQ[idx]) idx = i;
 
   switch(idx) {
     case 0:

src/FDM/JSBSim/math/FGMatrix33.h

@@ -49,7 +49,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_MATRIX33 "$Id: FGMatrix33.h,v 1.14 2010/12/07 12:57:14 jberndt Exp $"
+#define ID_MATRIX33 "$Id: FGMatrix33.h,v 1.15 2011/12/10 15:49:21 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -277,7 +277,7 @@ public:
   /** Returns the quaternion associated with this direction cosine (rotation) matrix.
   */
   FGQuaternion GetQuaternion(void);
-  
+
   /** Determinant of the matrix.
       @return the determinant of the matrix.
    */
@@ -465,7 +465,4 @@ std::ostream& operator<<(std::ostream& os, const FGMatrix33& M);
 std::istream& operator>>(std::istream& is, FGMatrix33& M);
 
 } // namespace JSBSim
-
-#include "FGQuaternion.h"
-
 #endif

src/FDM/JSBSim/math/FGModelFunctions.cpp

@@ -37,14 +37,16 @@ COMMENTS, REFERENCES,  and NOTES
 INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#include "FGModelFunctions.h"
+#include <iostream>
+#include <sstream>
 #include <string>
+#include "FGModelFunctions.h"
 
 using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGModelFunctions.cpp,v 1.4 2010/09/07 00:40:03 jberndt Exp $";
+static const char *IdSrc = "$Id: FGModelFunctions.cpp,v 1.5 2012/04/13 13:25:52 jberndt Exp $";
 static const char *IdHdr = ID_MODELFUNCTIONS;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -163,4 +165,54 @@ void FGModelFunctions::RunPostFunctions(void)
     (*it)->GetValue();
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+string FGModelFunctions::GetFunctionStrings(const string& delimeter) const
+{
+  string FunctionStrings = "";
+  bool firstime = true;
+  unsigned int sd;
+
+  for (sd = 0; sd < PreFunctions.size(); sd++) {
+    if (firstime) {
+      firstime = false;
+    } else {
+      FunctionStrings += delimeter;
+    }
+    FunctionStrings += PreFunctions[sd]->GetName();
+  }
+
+  for (sd = 0; sd < PostFunctions.size(); sd++) {
+    if (firstime) {
+      firstime = false;
+    } else {
+      FunctionStrings += delimeter;
+    }
+    FunctionStrings += PostFunctions[sd]->GetName();
+  }
+
+  return FunctionStrings;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+string FGModelFunctions::GetFunctionValues(const string& delimeter) const
+{
+  ostringstream buf;
+
+  for (unsigned int sd = 0; sd < PreFunctions.size(); sd++) {
+    if (buf.tellp() > 0) buf << delimeter;
+    buf << PreFunctions[sd]->GetValue();
+  }
+
+  for (unsigned int sd = 0; sd < PostFunctions.size(); sd++) {
+    if (buf.tellp() > 0) buf << delimeter;
+    buf << PostFunctions[sd]->GetValue();
+  }
+
+  return buf.str();
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 }

src/FDM/JSBSim/math/FGModelFunctions.h

@@ -44,7 +44,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_MODELFUNCTIONS "$Id: FGModelFunctions.h,v 1.4 2011/06/21 04:41:54 jberndt Exp $"
+#define ID_MODELFUNCTIONS "$Id: FGModelFunctions.h,v 1.5 2012/04/13 13:25:52 jberndt Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -83,6 +83,17 @@ public:
   void PreLoad(Element* el, FGPropertyManager* PropertyManager, std::string prefix="");
   void PostLoad(Element* el, FGPropertyManager* PropertyManager, std::string prefix="");
 
+  /** Gets the strings for the current set of functions.
+      @param delimeter either a tab or comma string depending on output type
+      @return a string containing the descriptive names for all functions */
+  std::string GetFunctionStrings(const std::string& delimeter) const;
+
+  /** Gets the function values.
+      @param delimeter either a tab or comma string depending on output type
+      @return a string containing the numeric values for the current set of
+      functions */
+  std::string GetFunctionValues(const std::string& delimeter) const;
+
 protected:
   std::vector <FGFunction*> PreFunctions;
   std::vector <FGFunction*> PostFunctions;

src/FDM/JSBSim/math/FGQuaternion.h

@@ -47,7 +47,7 @@ SENTRY
   DEFINITIONS
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_QUATERNION "$Id: FGQuaternion.h,v 1.23 2011/10/31 14:54:40 bcoconni Exp $"
+#define ID_QUATERNION "$Id: FGQuaternion.h,v 1.24 2011/12/10 15:49:21 bcoconni Exp $"
 
 namespace JSBSim {
 
@@ -463,6 +463,8 @@ public:
       Useful for initialization of increments */
   static FGQuaternion zero(void) { return FGQuaternion( 0.0, 0.0, 0.0, 0.0 ); }
 
+  friend FGQuaternion QExp(const FGColumnVector3& omega);
+
 private:
   /** Copying by assigning the vector valued components.  */
   FGQuaternion(double q1, double q2, double q3, double q4) : mCacheValid(false)
@@ -520,6 +522,24 @@ inline FGQuaternion operator*(double scalar, const FGQuaternion& q) {
   return FGQuaternion(scalar*q.data[0], scalar*q.data[1], scalar*q.data[2], scalar*q.data[3]);
 }
 
+/** Quaternion exponential
+    @param omega rotation velocity
+    Calculate the unit quaternion which is the result of the exponentiation of
+    the vector 'omega'.
+*/
+inline FGQuaternion QExp(const FGColumnVector3& omega) {
+  FGQuaternion qexp;
+  double angle = omega.Magnitude();
+  double sina_a = angle > 0.0 ? sin(angle)/angle : 1.0;
+
+  qexp.data[0] = cos(angle);
+  qexp.data[1] = omega(1) * sina_a;
+  qexp.data[2] = omega(2) * sina_a;
+  qexp.data[3] = omega(3) * sina_a;
+
+  return qexp;
+}
+
 /** Write quaternion to a stream.
     @param os Stream to write to.
     @param q Quaternion to write.

src/FDM/JSBSim/models/FGAccelerations.cpp

@@ -29,7 +29,7 @@
 FUNCTIONAL DESCRIPTION
 --------------------------------------------------------------------------------
 This class encapsulates the calculation of the derivatives of the state vectors
-UVW and PQR - the translational and rotational rates relative to the planet 
+UVW and PQR - the translational and rotational rates relative to the planet
 fixed frame. The derivatives relative to the inertial frame are also calculated
 as a side effect. Also, the derivative of the attitude quaterion is also calculated.
 
@@ -45,6 +45,8 @@ COMMENTS, REFERENCES,  and NOTES
 [2] Richard E. McFarland, "A Standard Kinematic Model for Flight Simulation at
     NASA-Ames", NASA CR-2497, January 1975
 [3] Erin Catto, "Iterative Dynamics with Temporal Coherence", February 22, 2005
+[4] Mark Harris and Robert Lyle, "Spacecraft Gravitational Torques",
+    NASA SP-8024, May 1969
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 INCLUDES
@@ -58,7 +60,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGAccelerations.cpp,v 1.8 2011/08/30 20:49:04 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGAccelerations.cpp,v 1.13 2012/02/18 19:11:37 bcoconni Exp $";
 static const char *IdHdr = ID_ACCELERATIONS;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -71,7 +73,8 @@ FGAccelerations::FGAccelerations(FGFDMExec* fdmex)
   Debug(0);
   Name = "FGAccelerations";
   gravType = gtWGS84;
- 
+  gravTorque = false;
+
   vPQRidot.InitMatrix();
   vUVWidot.InitMatrix();
   vGravAccel.InitMatrix();
@@ -112,16 +115,12 @@ bool FGAccelerations::Run(bool Holding)
   if (FGModel::Run(Holding)) return true;  // Fast return if we have nothing to do ...
   if (Holding) return false;
 
-  RunPreFunctions();
-
   CalculatePQRdot();   // Angular rate derivative
   CalculateUVWdot();   // Translational rate derivative
   CalculateQuatdot();  // Angular orientation derivative
 
   ResolveFrictionForces(in.DeltaT * rate);  // Update rate derivatives with friction forces
 
-  RunPostFunctions();
-
   Debug(2);
   return false;
 }
@@ -129,19 +128,29 @@ bool FGAccelerations::Run(bool Holding)
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Compute body frame rotational accelerations based on the current body moments
 //
-// vPQRdot is the derivative of the absolute angular velocity of the vehicle 
+// vPQRdot is the derivative of the absolute angular velocity of the vehicle
-// (body rate with respect to the inertial frame), expressed in the body frame,
+// (body rate with respect to the ECEF frame), expressed in the body frame,
 // where the derivative is taken in the body frame.
 // J is the inertia matrix
 // Jinv is the inverse inertia matrix
 // vMoments is the moment vector in the body frame
 // in.vPQRi is the total inertial angular velocity of the vehicle
 // expressed in the body frame.
-// Reference: See Stevens and Lewis, "Aircraft Control and Simulation", 
+// Reference: See Stevens and Lewis, "Aircraft Control and Simulation",
 //            Second edition (2004), eqn 1.5-16e (page 50)
 
 void FGAccelerations::CalculatePQRdot(void)
 {
+  if (gravTorque) {
+    // Compute the gravitational torque
+    // Reference: See Harris and Lyle "Spacecraft Gravitational Torques",
+    //            NASA SP-8024 (1969) eqn (2) (page 7)
+    FGColumnVector3 R = in.Ti2b * in.vInertialPosition;
+    double invRadius = 1.0 / R.Magnitude();
+    R *= invRadius;
+    in.Moment += (3.0 * in.GAccel * invRadius) * (R * (in.J * R));
+  }
+
   // Compute body frame rotational accelerations based on the current body
   // moments and the total inertial angular velocity expressed in the body
   // frame.
@@ -154,7 +163,7 @@ void FGAccelerations::CalculatePQRdot(void)
 // Compute the quaternion orientation derivative
 //
 // vQtrndot is the quaternion derivative.
-// Reference: See Stevens and Lewis, "Aircraft Control and Simulation", 
+// Reference: See Stevens and Lewis, "Aircraft Control and Simulation",
 //            Second edition (2004), eqn 1.5-16b (page 50)
 
 void FGAccelerations::CalculateQuatdot(void)
@@ -167,7 +176,7 @@ void FGAccelerations::CalculateQuatdot(void)
 // This set of calculations results in the body and inertial frame accelerations
 // being computed.
 // Compute body and inertial frames accelerations based on the current body
-// forces including centripetal and coriolis accelerations for the former.
+// forces including centripetal and Coriolis accelerations for the former.
 // in.vOmegaPlanet is the Earth angular rate - expressed in the inertial frame -
 //   so it has to be transformed to the body frame. More completely,
 //   in.vOmegaPlanet is the rate of the ECEF frame relative to the Inertial
@@ -177,7 +186,7 @@ void FGAccelerations::CalculateQuatdot(void)
 //   in the body frame.
 // in.vUVW is the vehicle velocity relative to the ECEF frame, expressed
 //   in the body frame.
-// Reference: See Stevens and Lewis, "Aircraft Control and Simulation", 
+// Reference: See Stevens and Lewis, "Aircraft Control and Simulation",
 //            Second edition (2004), eqns 1.5-13 (pg 48) and 1.5-16d (page 50)
 
 void FGAccelerations::CalculateUVWdot(void)
@@ -192,27 +201,30 @@ void FGAccelerations::CalculateUVWdot(void)
   // Include Gravitation accel
   switch (gravType) {
     case gtStandard:
-      vGravAccel = in.Tl2b * FGColumnVector3( 0.0, 0.0, in.GAccel );
+      {
+        double radius = in.vInertialPosition.Magnitude();
+        vGravAccel = -(in.GAccel / radius) * in.vInertialPosition;
+      }
       break;
     case gtWGS84:
-      vGravAccel = in.Tec2b * in.J2Grav;
+      vGravAccel = in.Tec2i * in.J2Grav;
       break;
   }
 
-  vUVWdot += vGravAccel;
+  vUVWdot += in.Ti2b * vGravAccel;
-  vUVWidot = in.Tb2i * (vBodyAccel + vGravAccel);
+  vUVWidot = in.Tb2i * vBodyAccel + vGravAccel;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Resolves the contact forces just before integrating the EOM.
 // This routine is using Lagrange multipliers and the projected Gauss-Seidel
 // (PGS) method.
-// Reference: See Erin Catto, "Iterative Dynamics with Temporal Coherence", 
+// Reference: See Erin Catto, "Iterative Dynamics with Temporal Coherence",
 //            February 22, 2005
 // In JSBSim there is only one rigid body (the aircraft) and there can be
 // multiple points of contact between the aircraft and the ground. As a
-// consequence our matrix J*M^-1*J^T is not sparse and the algorithm described
+// consequence our matrix Jac*M^-1*Jac^T is not sparse and the algorithm
-// in Catto's paper has been adapted accordingly.
+// described in Catto's paper has been adapted accordingly.
 // The friction forces are resolved in the body frame relative to the origin
 // (Earth center).
 
@@ -230,7 +242,7 @@ void FGAccelerations::ResolveFrictionForces(double dt)
   // If no gears are in contact with the ground then return
   if (!n) return;
 
-  vector<double> a(n*n); // Will contain J*M^-1*J^T
+  vector<double> a(n*n); // Will contain Jac*M^-1*Jac^T
   vector<double> rhs(n);
 
   // Assemble the linear system of equations
@@ -239,7 +251,7 @@ void FGAccelerations::ResolveFrictionForces(double dt)
     FGColumnVector3 v2 = Jinv * multipliers[i]->MomentJacobian; // Should be J^-T but J is symmetric and so is J^-1
 
     for (int j=0; j < i; j++)
-      a[i*n+j] = a[j*n+i]; // Takes advantage of the symmetry of J^T*M^-1*J
+      a[i*n+j] = a[j*n+i]; // Takes advantage of the symmetry of Jac^T*M^-1*Jac
     for (int j=i; j < n; j++)
       a[i*n+j] = DotProduct(v1, multipliers[j]->ForceJacobian)
                + DotProduct(v2, multipliers[j]->MomentJacobian);
@@ -249,12 +261,12 @@ void FGAccelerations::ResolveFrictionForces(double dt)
 
   // Translation
   vdot = vUVWdot;
-  if (dt > 0.) // Zeroes out the relative movement between aircraft and ground
+  if (dt > 0.) // Zeroes out the relative movement between the aircraft and the ground
     vdot += (in.vUVW - in.Tec2b * in.TerrainVelocity) / dt;
 
   // Rotation
   wdot = vPQRdot;
-  if (dt > 0.) // Zeroes out the relative movement between aircraft and ground
+  if (dt > 0.) // Zeroes out the relative movement between the aircraft and the ground
     wdot += (in.vPQR - in.Tec2b * in.TerrainAngularVel) / dt;
 
   // Prepare the linear system for the Gauss-Seidel algorithm :
@@ -277,7 +289,7 @@ void FGAccelerations::ResolveFrictionForces(double dt)
     for (int i=0; i < n; i++) {
       double lambda0 = multipliers[i]->value;
       double dlambda = rhs[i];
-      
+
       for (int j=0; j < n; j++)
         dlambda -= a[i*n+j]*multipliers[j]->value;
 
@@ -307,6 +319,7 @@ void FGAccelerations::ResolveFrictionForces(double dt)
   vPQRdot += omegadot;
   vPQRidot += omegadot;
 }
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 void FGAccelerations::InitializeDerivatives(void)
@@ -333,6 +346,7 @@ void FGAccelerations::bind(void)
   PropertyManager->Tie("accelerations/wdot-ft_sec2", this, eW, (PMF)&FGAccelerations::GetUVWdot);
 
   PropertyManager->Tie("simulation/gravity-model", &gravType);
+  PropertyManager->Tie("simulation/gravitational-torque", &gravTorque);
 
   PropertyManager->Tie("forces/fbx-total-lbs", this, eX, (PMF)&FGAccelerations::GetForces);
   PropertyManager->Tie("forces/fby-total-lbs", this, eY, (PMF)&FGAccelerations::GetForces);

src/FDM/JSBSim/models/FGAccelerations.h

@@ -50,7 +50,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_ACCELERATIONS "$Id: FGAccelerations.h,v 1.8 2011/10/31 14:54:41 bcoconni Exp $"
+#define ID_ACCELERATIONS "$Id: FGAccelerations.h,v 1.12 2012/02/19 14:07:27 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -64,13 +64,40 @@ CLASS DOCUMENTATION
 
 /** Handles the calculation of accelerations.
 
-    -Calculate the angular accelerations
+    - Calculate the angular accelerations
-    -Calculate the translational accelerations
+    - Calculate the translational accelerations
-    -Calculate the angular rate
+    - Calculate the angular rate
-    -Calculate the translational velocity
+    - Calculate the translational velocity
+
+    This class is collecting all the forces and the moments acting on the body
+    to calculate the corresponding accelerations according to Newton's second
+    law. This is also where the friction forces related to the ground reactions
+    are evaluated.
+
+    JSBSim provides several ways to calculate the influence of the gravity on
+    the vehicle. The different options can be selected via the following
+    properties :
+    @property simulation/gravity-model (read/write) Selects the gravity model.
+              Two options are available : 0 (Standard gravity assuming the Earth
+              is spherical) or 1 (WGS84 gravity taking the Earth oblateness into
+              account). WGS84 gravity is the default.
+    @property simulation/gravitational-torque (read/write) Enables/disables the
+              calculations of the gravitational torque on the vehicle. This is
+              mainly relevant for spacecrafts that are orbiting at low altitudes.
+              Gravitational torque calculations are disabled by default.
+
+    Special care is taken in the calculations to obtain maximum fidelity in
+    JSBSim results. In FGAccelerations, this is obtained by avoiding as much as
+    possible the transformations from one frame to another. As a consequence,
+    the frames in which the accelerations are primarily evaluated are dictated
+    by the frames in which FGPropagate resolves the equations of movement (the
+    ECI frame for the translations and the body frame for the rotations).
+
+    @see Mark Harris and Robert Lyle, "Spacecraft Gravitational Torques",
+         NASA SP-8024, May 1969
 
     @author Jon S. Berndt, Mathias Froehlich, Bertrand Coconnier
-    @version $Id: FGAccelerations.h,v 1.8 2011/10/31 14:54:41 bcoconni Exp $
+    @version $Id: FGAccelerations.h,v 1.12 2012/02/19 14:07:27 bcoconni Exp $
   */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -85,30 +112,44 @@ public:
 
   /// Destructor
   ~FGAccelerations();
-  
+
   /// These define the indices use to select the gravitation models.
-  enum eGravType {gtStandard, gtWGS84}; 
+  enum eGravType {
+    /// Evaluate gravity using Newton's classical formula assuming the Earth is spherical
+    gtStandard,
+    /// Evaluate gravity using WGS84 formulas that take the Earth oblateness into account
+    gtWGS84
+  };
 
   /** Initializes the FGAccelerations class after instantiation and prior to first execution.
-      The base class FGModel::InitModel is called first, initializing pointers to the 
+      The base class FGModel::InitModel is called first, initializing pointers to the
       other FGModel objects (and others).  */
   bool InitModel(void);
 
   /** Runs the state propagation model; called by the Executive
       Can pass in a value indicating if the executive is directing the simulation to Hold.
-      @param Holding if true, the executive has been directed to hold the sim from 
+      @param Holding if true, the executive has been directed to hold the sim from
                      advancing time. Some models may ignore this flag, such as the Input
                      model, which may need to be active to listen on a socket for the
                      "Resume" command to be given.
       @return false if no error */
   bool Run(bool Holding);
 
+  /** Retrieves the time derivative of the body orientation quaternion.
+      Retrieves the time derivative of the body orientation quaternion based on
+      the rate of change of the orientation between the body and the ECI frame.
+      The quaternion returned is represented by an FGQuaternion reference. The
+      quaternion is 1-based, so that the first element can be retrieved using
+      the "()" operator.
+      units rad/sec^2
+      @return The time derivative of the body orientation quaternion.
+  */
   const FGQuaternion& GetQuaterniondot(void) const {return vQtrndot;}
 
   /** Retrieves the body axis acceleration.
       Retrieves the computed body axis accelerations based on the
       applied forces and accounting for a rotating body frame.
-      The vector returned is represented by an FGColumnVector reference. The vector
+      The vector returned is represented by an FGColumnVector3 reference. The vector
       for the acceleration in Body frame is organized (Ax, Ay, Az). The vector
       is 1-based, so that the first element can be retrieved using the "()" operator.
       In other words, vUVWdot(1) is Ax. Various convenience enumerators are defined
@@ -119,13 +160,27 @@ public:
   */
   const FGColumnVector3& GetUVWdot(void) const { return vUVWdot; }
 
+  /** Retrieves the body axis acceleration in the ECI frame.
+      Retrieves the computed body axis accelerations based on the applied forces.
+      The ECI frame being an inertial frame this vector does not contain the
+      Coriolis and centripetal accelerations. The vector is expressed in the
+      Body frame.
+      The vector returned is represented by an FGColumnVector3 reference. The
+      vector for the acceleration in Body frame is organized (Aix, Aiy, Aiz). The
+      vector is 1-based, so that the first element can be retrieved using the
+      "()" operator. In other words, vUVWidot(1) is Aix. Various convenience
+      enumerators are defined in FGJSBBase. The relevant enumerators for the
+      vector returned by this call are, eX=1, eY=2, eZ=3.
+      units ft/sec^2
+      @return Body axis translational acceleration in ft/sec^2.
+  */
   const FGColumnVector3& GetUVWidot(void) const { return vUVWidot; }
 
   /** Retrieves the body axis angular acceleration vector.
       Retrieves the body axis angular acceleration vector in rad/sec^2. The
-      angular acceleration vector is determined from the applied forces and
+      angular acceleration vector is determined from the applied moments and
       accounts for a rotating frame.
-      The vector returned is represented by an FGColumnVector reference. The vector
+      The vector returned is represented by an FGColumnVector3 reference. The vector
       for the angular acceleration in Body frame is organized (Pdot, Qdot, Rdot). The vector
       is 1-based, so that the first element can be retrieved using the "()" operator.
       In other words, vPQRdot(1) is Pdot. Various convenience enumerators are defined
@@ -135,12 +190,25 @@ public:
       @return The angular acceleration vector.
   */
   const FGColumnVector3& GetPQRdot(void) const {return vPQRdot;}
-  
+
+  /** Retrieves the axis angular acceleration vector in the ECI frame.
+      Retrieves the body axis angular acceleration vector measured in the ECI
+      frame and expressed in the body frame. The angular acceleration vector is
+      determined from the applied moments.
+      The vector returned is represented by an FGColumnVector3 reference. The
+      vector for the angular acceleration in Body frame is organized (Pidot,
+      Qidot, Ridot). The vector is 1-based, so that the first element can be
+      retrieved using the "()" operator. In other words, vPQRidot(1) is Pidot.
+      Various convenience enumerators are defined in FGJSBBase. The relevant
+      enumerators for the vector returned by this call are, eP=1, eQ=2, eR=3.
+      units rad/sec^2
+      @return The angular acceleration vector.
+  */
   const FGColumnVector3& GetPQRidot(void) const {return vPQRidot;}
 
   /** Retrieves a body frame acceleration component.
       Retrieves a body frame acceleration component. The acceleration returned
-      is extracted from the vUVWdot vector (an FGColumnVector). The vector for
+      is extracted from the vUVWdot vector (an FGColumnVector3). The vector for
       the acceleration in Body frame is organized (Ax, Ay, Az). The vector is
       1-based. In other words, GetUVWdot(1) returns Ax. Various convenience
       enumerators are defined in FGJSBBase. The relevant enumerators for the
@@ -151,14 +219,39 @@ public:
   */
   double GetUVWdot(int idx) const { return vUVWdot(idx); }
 
+  /** Retrieves the acceleration resulting from the applied forces.
+      Retrieves the ratio of the sum of all forces applied on the craft to its
+      mass. This does include the friction forces but not the gravity.
+      The vector returned is represented by an FGColumnVector3 reference. The
+      vector for the acceleration in Body frame is organized (Ax, Ay, Az). The
+      vector is 1-based, so that the first element can be retrieved using the
+      "()" operator. In other words, vBodyAccel(1) is Ax. Various convenience
+      enumerators are defined in FGJSBBase. The relevant enumerators for the
+      vector returned by this call are, eX=1, eY=2, eZ=3.
+      units ft/sec^2
+      @return The acceleration resulting from the applied forces.
+  */
   const FGColumnVector3& GetBodyAccel(void) const { return vBodyAccel; }
 
+
+  /** Retrieves a component of the acceleration resulting from the applied forces.
+      Retrieves a component of the ratio between the sum of all forces applied
+      on the craft to its mass. The value returned is extracted from the vBodyAccel
+      vector (an FGColumnVector3). The vector for the acceleration in Body frame
+      is organized (Ax, Ay, Az). The vector is 1-based. In other words,
+      GetBodyAccel(1) returns Ax. Various convenience enumerators are defined
+      in FGJSBBase. The relevant enumerators for the vector returned by this
+      call are, eX=1, eY=2, eZ=3.
+      units ft/sec^2
+      @param idx the index of the acceleration component desired (1-based).
+      @return The component of the acceleration resulting from the applied forces.
+  */
   double GetBodyAccel(int idx) const { return vBodyAccel(idx); }
 
   /** Retrieves a body frame angular acceleration component.
       Retrieves a body frame angular acceleration component. The angular
       acceleration returned is extracted from the vPQRdot vector (an
-      FGColumnVector). The vector for the angular acceleration in Body frame
+      FGColumnVector3). The vector for the angular acceleration in Body frame
       is organized (Pdot, Qdot, Rdot). The vector is 1-based. In other words,
       GetPQRdot(1) returns Pdot (roll acceleration). Various convenience
       enumerators are defined in FGJSBBase. The relevant enumerators for the
@@ -169,38 +262,111 @@ public:
   */
   double GetPQRdot(int axis) const {return vPQRdot(axis);}
 
+  /** Retrieves a component of the total moments applied on the body.
+      Retrieves a component of the total moments applied on the body. This does
+      include the moments generated by friction forces and the gravitational
+      torque (if the property \e simulation/gravitational-torque is set to true).
+      The vector for the total moments in the body frame is organized (Mx, My
+      , Mz). The vector is 1-based. In other words, GetMoments(1) returns Mx.
+      Various convenience enumerators are defined in FGJSBBase. The relevant
+      enumerators for the moments returned by this call are, eX=1, eY=2, eZ=3.
+      units lbs*ft
+      @param idx the index of the moments component desired (1-based).
+      @return The total moments applied on the body.
+   */
   double GetMoments(int idx) const { return in.Moment(idx) + vFrictionMoments(idx); }
+
+  /** Retrieves the total forces applied on the body.
+      Retrieves the total forces applied on the body. This does include the
+      friction forces but not the gravity.
+      The vector for the total forces in the body frame is organized (Fx, Fy
+      , Fz). The vector is 1-based. In other words, GetForces(1) returns Fx.
+      Various convenience enumerators are defined in FGJSBBase. The relevant
+      enumerators for the forces returned by this call are, eX=1, eY=2, eZ=3.
+      units lbs
+      @param idx the index of the forces component desired (1-based).
+      @return The total forces applied on the body.
+   */
   double GetForces(int idx) const { return in.Force(idx) + vFrictionForces(idx); }
+
+  /** Retrieves the ground moments applied on the body.
+      Retrieves the ground moments applied on the body. This does include the
+      ground normal reaction and friction moments.
+      The vector for the ground moments in the body frame is organized (Mx, My
+      , Mz). The vector is 1-based. In other words, GetGroundMoments(1) returns
+      Mx. Various convenience enumerators are defined in FGJSBBase. The relevant
+      enumerators for the moments returned by this call are, eX=1, eY=2, eZ=3.
+      units lbs*ft
+      @param idx the index of the moments component desired (1-based).
+      @return The ground moments applied on the body.
+   */
   double GetGroundMoments(int idx) const { return in.GroundMoment(idx) + vFrictionMoments(idx); }
+
+  /** Retrieves the ground forces applied on the body.
+      Retrieves the ground forces applied on the body. This does include the
+      ground normal reaction and friction forces.
+      The vector for the total moments in the body frame is organized (Fx, Fy
+      , Fz). The vector is 1-based. In other words, GetGroundForces(1) returns
+      Fx. Various convenience enumerators are defined in FGJSBBase. The relevant
+      enumerators for the forces returned by this call are, eX=1, eY=2, eZ=3.
+      units lbs
+      @param idx the index of the forces component desired (1-based).
+      @return The ground forces applied on the body.
+   */
   double GetGroundForces(int idx) const { return in.GroundForce(idx) + vFrictionForces(idx); }
 
+  /** Initializes the FGAccelerations class prior to a new execution.
+      Initializes the class prior to a new execution when the input data stored
+      in the Inputs structure have been set to their initial values.
+   */
   void InitializeDerivatives(void);
 
-  void DumpState(void);
-
   struct Inputs {
+    /// The body inertia matrix expressed in the body frame
     FGMatrix33 J;
+    /// The inverse of the inertia matrix J
     FGMatrix33 Jinv;
+    /// Transformation matrix from the ECI to the Body frame
     FGMatrix33 Ti2b;
+    /// Transformation matrix from the Body to the ECI frame
     FGMatrix33 Tb2i;
+    /// Transformation matrix from the ECEF to the Body frame
     FGMatrix33 Tec2b;
-    FGMatrix33 Tl2b;
+    /// Transformation matrix from the ECEF to the ECI frame
+    FGMatrix33 Tec2i;
+    /// Orientation quaternion of the body with respect to the ECI frame
     FGQuaternion qAttitudeECI;
+    /// Total moments applied to the body except friction and gravity (expressed in the body frame)
     FGColumnVector3 Moment;
+    /// Moments generated by the ground normal reactions expressed in the body frame. Does not account for friction.
     FGColumnVector3 GroundMoment;
+    /// Total forces applied to the body except friction and gravity (expressed in the body frame)
     FGColumnVector3 Force;
+    /// Forces generated by the ground normal reactions expressed in the body frame. Does not account for friction.
     FGColumnVector3 GroundForce;
+    /// Gravity intensity vector using WGS84 formulas (expressed in the ECEF frame).
     FGColumnVector3 J2Grav;
+    /// Angular velocities of the body with respect to the ECI frame (expressed in the body frame).
     FGColumnVector3 vPQRi;
+    /// Angular velocities of the body with respect to the local frame (expressed in the body frame).
     FGColumnVector3 vPQR;
+    /// Velocities of the body with respect to the local frame (expressed in the body frame).
     FGColumnVector3 vUVW;
+    /// Body position (X,Y,Z) measured in the ECI frame.
     FGColumnVector3 vInertialPosition;
+    /// Earth rotating vector (expressed in the ECI frame).
     FGColumnVector3 vOmegaPlanet;
+    /// Terrain velocities with respect to the local frame (expressed in the ECEF frame).
     FGColumnVector3 TerrainVelocity;
+    /// Terrain angular velocities with respect to the local frame (expressed in the ECEF frame).
     FGColumnVector3 TerrainAngularVel;
+    /// Time step
     double DeltaT;
+    /// Body mass
     double Mass;
+    /// Gravity intensity assuming the Earth is spherical
     double GAccel;
+    /// List of Lagrange multipliers set by FGLGear for friction forces calculations.
     std::vector<LagrangeMultiplier*> *MultipliersList;
   } in;
 
@@ -215,6 +381,7 @@ private:
   FGColumnVector3 vFrictionMoments;
 
   int gravType;
+  bool gravTorque;
 
   void CalculatePQRdot(void);
   void CalculateQuatdot(void);

src/FDM/JSBSim/models/FGAerodynamics.cpp

@@ -48,7 +48,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGAerodynamics.cpp,v 1.44 2011/10/23 14:23:13 jentron Exp $";
+static const char *IdSrc = "$Id: FGAerodynamics.cpp,v 1.45 2012/04/13 13:25:52 jberndt Exp $";
 static const char *IdHdr = ID_AERODYNAMICS;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -372,6 +372,17 @@ string FGAerodynamics::GetAeroFunctionStrings(const string& delimeter) const
       AeroFunctionStrings += AeroFunctions[axis][sd]->GetName();
     }
   }
+
+  string FunctionStrings = FGModelFunctions::GetFunctionStrings(delimeter);
+
+  if (FunctionStrings.size() > 0) {
+    if (AeroFunctionStrings.size() > 0) {
+      AeroFunctionStrings += delimeter + FunctionStrings;
+    } else {
+      AeroFunctionStrings = FunctionStrings;
+    }
+  }
+
   return AeroFunctionStrings;
 }
 
@@ -388,6 +399,16 @@ string FGAerodynamics::GetAeroFunctionValues(const string& delimeter) const
     }
   }
 
+  string FunctionValues = FGModelFunctions::GetFunctionValues(delimeter);
+
+  if (FunctionValues.size() > 0) {
+    if (buf.str().size() > 0) {
+      buf << delimeter << FunctionValues;
+    } else {
+      buf << FunctionValues;
+    }
+  }
+
   return buf.str();
 }
 

src/FDM/JSBSim/models/FGAtmosphere.cpp

@@ -50,7 +50,7 @@ INCLUDES
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGAtmosphere.cpp,v 1.49 2011/09/11 11:36:04 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGAtmosphere.cpp,v 1.51 2012/04/13 13:18:28 jberndt Exp $";
 static const char *IdHdr = ID_ATMOSPHERE;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -101,12 +101,8 @@ bool FGAtmosphere::Run(bool Holding)
   if (FGModel::Run(Holding)) return true;
   if (Holding) return false;
 
-  RunPreFunctions();
-
   Calculate(in.altitudeASL);
 
-  RunPostFunctions();
-
   Debug(2);
   return false;
 }
@@ -128,7 +124,7 @@ void FGAtmosphere::Calculate(double altitude)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-void FGAtmosphere::SetPressureSL(double pressure, ePressure unit)
+void FGAtmosphere::SetPressureSL(ePressure unit, double pressure)
 {
   double press = ConvertToPSF(pressure, unit);
 
@@ -217,7 +213,9 @@ void FGAtmosphere::bind(void)
   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/P-sl-psf", this, ePSF,
+//                                   (PMFi)&FGAtmosphere::GetPressureSL,
+//                                   (PMF)&FGAtmosphere::SetPressureSL);
   PropertyManager->Tie("atmosphere/a-sl-fps", this, &FGAtmosphere::GetSoundSpeedSL);
   PropertyManager->Tie("atmosphere/theta", this, &FGAtmosphere::GetTemperatureRatio);
   PropertyManager->Tie("atmosphere/sigma", this, &FGAtmosphere::GetDensityRatio);
@@ -264,7 +262,7 @@ void FGAtmosphere::Debug(int from)
   }
   if (debug_lvl & 16) { // Sanity checking
   }
-  if (debug_lvl & 128) { // 
+  if (debug_lvl & 128) { //
   }
   if (debug_lvl & 64) {
     if (from == 0) { // Constructor

src/FDM/JSBSim/models/FGAtmosphere.h

@@ -45,7 +45,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_ATMOSPHERE "$Id: FGAtmosphere.h,v 1.29 2011/07/10 20:18:14 jberndt Exp $"
+#define ID_ATMOSPHERE "$Id: FGAtmosphere.h,v 1.30 2012/04/13 13:18:28 jberndt Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -74,7 +74,7 @@ CLASS DOCUMENTATION
   @property atmosphere/a-ratio
 
   @author Jon Berndt
-  @version $Id: FGAtmosphere.h,v 1.29 2011/07/10 20:18:14 jberndt Exp $
+  @version $Id: FGAtmosphere.h,v 1.30 2012/04/13 13:18:28 jberndt Exp $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -161,10 +161,10 @@ public:
   virtual double GetPressureRatio(void) const { return Pressure*rSLpressure; }
 
   /** Sets the sea level pressure for modeling.
-      @param pressure The pressure in the units specified (PSF by default).
+      @param pressure The pressure in the units specified.
       @param unit the unit of measure that the specified pressure is
                   supplied in.*/
-  virtual void SetPressureSL(double pressure, ePressure unit=ePSF);
+  virtual void SetPressureSL(ePressure unit, double pressure);
   //@}
 
   //  *************************************************************************

src/FDM/JSBSim/models/FGAuxiliary.cpp

@@ -50,7 +50,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGAuxiliary.cpp,v 1.55 2011/11/12 18:59:11 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGAuxiliary.cpp,v 1.56 2011/12/11 17:03:05 bcoconni Exp $";
 static const char *IdHdr = ID_AUXILIARY;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -124,7 +124,7 @@ bool FGAuxiliary::InitModel(void)
 
   return true;
 }
-  
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 FGAuxiliary::~FGAuxiliary()
@@ -141,9 +141,7 @@ bool FGAuxiliary::Run(bool Holding)
   if (FGModel::Run(Holding)) return true; // return true if error returned from base class
   if (Holding) return false;
 
-  RunPreFunctions();
+  // Rotation
-
-// Rotation
 
   vEulerRates(eTht) = in.vPQR(eQ)*in.CosPhi - in.vPQR(eR)*in.SinPhi;
   if (in.CosTht != 0.0) {
@@ -151,7 +149,7 @@ bool FGAuxiliary::Run(bool Holding)
     vEulerRates(ePhi) = in.vPQR(eP) + vEulerRates(ePsi)*in.SinTht;
   }
 
-// Combine the wind speed with aircraft speed to obtain wind relative speed
+  // Combine the wind speed with aircraft speed to obtain wind relative speed
   vAeroPQR = in.vPQR - in.TurbPQR;
   vAeroUVW = in.vUVW - in.Tl2b * in.TotalWindNED;
 
@@ -195,7 +193,7 @@ bool FGAuxiliary::Run(bool Holding)
   vMachUVW(eW) = vAeroUVW(eW) / in.SoundSpeed;
   double MachU2 = MachU * MachU;
 
-// Position
+  // Position
 
   Vground = sqrt( in.vVel(eNorth)*in.vVel(eNorth) + in.vVel(eEast)*in.vVel(eEast) );
 
@@ -259,8 +257,6 @@ bool FGAuxiliary::Run(bool Holding)
   // When all models are executed calculate the distance from the initial point.
   CalculateRelativePosition();
 
-  RunPostFunctions();
-
   return false;
 }
 
@@ -344,7 +340,7 @@ double FGAuxiliary::GetNlf(void) const
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 void FGAuxiliary::CalculateRelativePosition(void)  //ToDo: This belongs elsewhere - perhaps in FGPropagate or Exec
-{ 
+{
   const double earth_radius_mt = in.ReferenceRadius*fttom;
   lat_relative_position=(in.Latitude  - FDMExec->GetIC()->GetLatitudeDegIC() *degtorad)*earth_radius_mt;
   lon_relative_position=(in.Longitude - FDMExec->GetIC()->GetLongitudeDegIC()*degtorad)*earth_radius_mt;

src/FDM/JSBSim/models/FGGroundReactions.cpp

@@ -46,7 +46,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGGroundReactions.cpp,v 1.36 2011/08/21 15:13:22 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGGroundReactions.cpp,v 1.39 2012/04/01 17:05:51 bcoconni Exp $";
 static const char *IdHdr = ID_GROUNDREACTIONS;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -140,8 +140,6 @@ bool FGGroundReactions::Load(Element* el)
 
   FGModel::Load(el); // Perform base class Load
 
-  in.vWhlBodyVec.resize(lGear.size());
-
   for (unsigned int i=0; i<lGear.size();i++) lGear[i]->bind();
 
   PostLoad(el, PropertyManager);

src/FDM/JSBSim/models/FGInertial.cpp

@@ -43,7 +43,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGInertial.cpp,v 1.25 2011/10/31 14:54:41 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGInertial.cpp,v 1.26 2011/12/11 17:03:05 bcoconni Exp $";
 static const char *IdHdr = ID_INERTIAL;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -61,7 +61,7 @@ FGInertial::FGInertial(FGFDMExec* fgex) : FGModel(fgex)
   RadiusReference = 20925650.00;        // Equatorial radius (WGS84)
   C2_0            = -4.84165371736E-04; // WGS84 value for the C2,0 coefficient
   J2              = 1.0826266836E-03;   // WGS84 value for J2
-  a               = 20925646.3255;      // WGS84 semimajor axis length in feet 
+  a               = 20925646.3255;      // WGS84 semimajor axis length in feet
   b               = 20855486.5951;      // WGS84 semiminor axis length in feet
 
   // Lunar defaults
@@ -71,7 +71,7 @@ FGInertial::FGInertial(FGFDMExec* fgex) : FGModel(fgex)
   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 
+  a               = 5702559.05;           // semimajor axis length in feet
   b               = 5695439.63;           // semiminor axis length in feet
   */
 
@@ -106,13 +106,9 @@ bool FGInertial::Run(bool Holding)
   if (FGModel::Run(Holding)) return true;
   if (Holding) return false;
 
-  RunPreFunctions();
-
   // Gravitation accel
   gAccel = GetGAccel(in.Radius);
 
-  RunPostFunctions();
-
   return false;
 }
 

src/FDM/JSBSim/models/FGInput.cpp

@@ -53,7 +53,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGInput.cpp,v 1.21 2011/05/20 03:18:36 jberndt Exp $";
+static const char *IdSrc = "$Id: FGInput.cpp,v 1.22 2012/01/21 16:46:09 jberndt Exp $";
 static const char *IdHdr = ID_INPUT;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -176,6 +176,22 @@ bool FGInput::Run(bool Holding)
 
         FDMExec->Resume();
         socket->Reply("");
+	
+      } else if (command == "iterate") {             // ITERATE
+
+        int argumentInt;
+        istringstream (argument) >> argumentInt;
+        if (argument.size() == 0) {
+          socket->Reply("No argument supplied for number of iterations.\n");
+          break;
+        }
+        if ( !(argumentInt > 0) ){
+          socket->Reply("Required argument must be a positive Integer.\n");
+          break;
+        }
+        FDMExec->EnableIncrementThenHold( argumentInt );
+        FDMExec->Resume();
+        socket->Reply("");
 
       } else if (command == "quit") {                   // QUIT
 

src/FDM/JSBSim/models/FGLGear.cpp

@@ -60,12 +60,13 @@ DEFINITIONS
 GLOBAL DATA
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-static const char *IdSrc = "$Id: FGLGear.cpp,v 1.92 2011/11/10 12:06:14 jberndt Exp $";
+static const char *IdSrc = "$Id: FGLGear.cpp,v 1.100 2012/04/01 17:05:51 bcoconni Exp $";
 static const char *IdHdr = ID_LGEAR;
 
 // Body To Structural (body frame is rotated 180 deg about Y and lengths are given in
 // ft instead of inches)
 const FGMatrix33 FGLGear::Tb2s(-1./inchtoft, 0., 0., 0., 1./inchtoft, 0., 0., 0., -1./inchtoft);
+const FGMatrix33 FGLGear::Ts2b(-inchtoft, 0., 0., 0., inchtoft, 0., 0., 0., -inchtoft);
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CLASS IMPLEMENTATION
@@ -79,19 +80,13 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
   Castered(false),
   StaticFriction(false)
 {
-  Element *force_table=0;
-  Element *dampCoeff=0;
-  Element *dampCoeffRebound=0;
-  string force_type="";
-
   kSpring = bDamp = bDampRebound = dynamicFCoeff = staticFCoeff = rollingFCoeff = maxSteerAngle = 0;
-  sSteerType = sBrakeGroup = sSteerType = "";
+  isRetractable = false;
-  isRetractable = 0;
   eDampType = dtLinear;
   eDampTypeRebound = dtLinear;
 
   name = el->GetAttributeValue("name");
-  sContactType = el->GetAttributeValue("type");
+  string sContactType = el->GetAttributeValue("type");
   if (sContactType == "BOGEY") {
     eContactType = ctBOGEY;
   } else if (sContactType == "STRUCTURE") {
@@ -101,7 +96,7 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
     eContactType = ctSTRUCTURE;
   }
 
-  // Default values for structural contact points 
+  // Default values for structural contact points
   if (eContactType == ctSTRUCTURE) {
     kSpring = in.EmptyWeight;
     bDamp = kSpring;
@@ -113,7 +108,7 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
   if (el->FindElement("spring_coeff"))
     kSpring = el->FindElementValueAsNumberConvertTo("spring_coeff", "LBS/FT");
   if (el->FindElement("damping_coeff")) {
-    dampCoeff = el->FindElement("damping_coeff");
+    Element* dampCoeff = el->FindElement("damping_coeff");
     if (dampCoeff->GetAttributeValue("type") == "SQUARE") {
       eDampType = dtSquare;
       bDamp   = el->FindElementValueAsNumberConvertTo("damping_coeff", "LBS/FT2/SEC2");
@@ -123,7 +118,7 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
   }
 
   if (el->FindElement("damping_coeff_rebound")) {
-    dampCoeffRebound = el->FindElement("damping_coeff_rebound");
+    Element* dampCoeffRebound = el->FindElement("damping_coeff_rebound");
     if (dampCoeffRebound->GetAttributeValue("type") == "SQUARE") {
       eDampTypeRebound = dtSquare;
       bDampRebound   = el->FindElementValueAsNumberConvertTo("damping_coeff_rebound", "LBS/FT2/SEC2");
@@ -141,32 +136,38 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
     staticFCoeff = el->FindElementValueAsNumber("static_friction");
   if (el->FindElement("rolling_friction"))
     rollingFCoeff = el->FindElementValueAsNumber("rolling_friction");
-  if (el->FindElement("max_steer"))
-    maxSteerAngle = el->FindElementValueAsNumberConvertTo("max_steer", "DEG");
   if (el->FindElement("retractable"))
     isRetractable = ((unsigned int)el->FindElementValueAsNumber("retractable"))>0.0?true:false;
 
+  if (el->FindElement("max_steer"))
+    maxSteerAngle = el->FindElementValueAsNumberConvertTo("max_steer", "DEG");
+
+  if (maxSteerAngle == 360) {
+    eSteerType = stCaster;
+    Castered = true;
+  }
+  else if (maxSteerAngle == 0.0) {
+    eSteerType = stFixed;
+  }
+  else
+    eSteerType = stSteer;
+
   GroundReactions = fdmex->GetGroundReactions();
   PropertyManager = fdmex->GetPropertyManager();
 
   ForceY_Table = 0;
-  force_table = el->FindElement("table");
+  Element* force_table = el->FindElement("table");
   while (force_table) {
-    force_type = force_table->GetAttributeValue("type");
+    string force_type = force_table->GetAttributeValue("type");
     if (force_type == "CORNERING_COEFF") {
       ForceY_Table = new FGTable(PropertyManager, force_table);
+      break;
     } else {
       cerr << "Undefined force table for " << name << " contact point" << endl;
     }
     force_table = el->FindNextElement("table");
   }
 
-  sBrakeGroup = el->FindElementValue("brake_group");
-
-  if (maxSteerAngle == 360) sSteerType = "CASTERED";
-  else if (maxSteerAngle == 0.0) sSteerType = "FIXED";
-  else sSteerType = "STEERABLE";
-
   Element* element = el->FindElement("location");
   if (element) vXYZn = element->FindElementTripletConvertTo("IN");
   else {cerr << "No location given for contact " << name << endl; exit(-1);}
@@ -174,25 +175,9 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
 
   element = el->FindElement("orientation");
   if (element && (eContactType == ctBOGEY)) {
-    vGearOrient = element->FindElementTripletConvertTo("RAD");
+    FGQuaternion quatFromEuler(element->FindElementTripletConvertTo("RAD"));
 
-    double cp,sp,cr,sr,cy,sy;
+    mTGear = quatFromEuler.GetT();
-
-    cp=cos(vGearOrient(ePitch)); sp=sin(vGearOrient(ePitch));
-    cr=cos(vGearOrient(eRoll));  sr=sin(vGearOrient(eRoll));
-    cy=cos(vGearOrient(eYaw));   sy=sin(vGearOrient(eYaw));
-
-    mTGear(1,1) =  cp*cy;
-    mTGear(2,1) =  cp*sy;
-    mTGear(3,1) = -sp;
-
-    mTGear(1,2) = sr*sp*cy - cr*sy;
-    mTGear(2,2) = sr*sp*sy + cr*cy;
-    mTGear(3,2) = sr*cp;
-
-    mTGear(1,3) = cr*sp*cy + sr*sy;
-    mTGear(2,3) = cr*sp*sy - sr*cy;
-    mTGear(3,3) = cr*cp;
   }
   else {
     mTGear(1,1) = 1.;
@@ -200,34 +185,22 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
     mTGear(3,3) = 1.;
   }
 
+  string sBrakeGroup = el->FindElementValue("brake_group");
+
   if      (sBrakeGroup == "LEFT"  ) eBrakeGrp = bgLeft;
   else if (sBrakeGroup == "RIGHT" ) eBrakeGrp = bgRight;
   else if (sBrakeGroup == "CENTER") eBrakeGrp = bgCenter;
-  else if (sBrakeGroup == "NOSE"  ) eBrakeGrp = bgNose;
+  else if (sBrakeGroup == "NOSE"  ) eBrakeGrp = bgCenter; // Nose brake is not supported by FGFCS
-  else if (sBrakeGroup == "TAIL"  ) eBrakeGrp = bgTail;
+  else if (sBrakeGroup == "TAIL"  ) eBrakeGrp = bgCenter; // Tail brake is not supported by FGFCS
   else if (sBrakeGroup == "NONE"  ) eBrakeGrp = bgNone;
-  else if (sBrakeGroup.empty()    ) {eBrakeGrp = bgNone;
+  else if (sBrakeGroup.empty()    ) eBrakeGrp = bgNone;
-                                     sBrakeGroup = "NONE (defaulted)";}
   else {
     cerr << "Improper braking group specification in config file: "
          << sBrakeGroup << " is undefined." << endl;
   }
 
-  if      (sSteerType == "STEERABLE") eSteerType = stSteer;
-  else if (sSteerType == "FIXED"    ) eSteerType = stFixed;
-  else if (sSteerType == "CASTERED" ) {eSteerType = stCaster; Castered = true;}
-  else if (sSteerType.empty()       ) {eSteerType = stFixed;
-                                       sSteerType = "FIXED (defaulted)";}
-  else {
-    cerr << "Improper steering type specification in config file: "
-         << sSteerType << " is undefined." << endl;
-  }
-
-  GearUp = false;
-  GearDown = true;
   GearPos  = 1.0;
   useFCSGearPos = false;
-  Servicable = true;
 
 // Add some AI here to determine if gear is located properly according to its
 // brake group type ??
@@ -245,11 +218,9 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs&
 
   compressLength  = 0.0;
   compressSpeed   = 0.0;
-  brakePct        = 0.0;
   maxCompLen      = 0.0;
 
   WheelSlip = 0.0;
-  TirePressureNorm = 1.0;
 
   // Set Pacejka terms
 
@@ -276,56 +247,51 @@ FGLGear::~FGLGear()
 
 const FGColumnVector3& FGLGear::GetBodyForces(void)
 {
+  double gearPos = 1.0;
   double t = fdmex->GetSimTime();
 
   vFn.InitMatrix();
 
-  if (isRetractable) ComputeRetractionState();
+  if (isRetractable) gearPos = GetGearUnitPos();
 
-  if (GearDown) {
+  if (gearPos > 0.99) { // Gear DOWN
-    FGColumnVector3 terrainVel, dummy;
+    FGColumnVector3 normal, terrainVel, dummy;
-
+    FGLocation gearLoc, contact;
-    vLocalGear = in.Tb2l * in.vWhlBodyVec[GearNumber]; // Get local frame wheel location
+    FGColumnVector3 vWhlBodyVec = Ts2b * (vXYZn - in.vXYZcg);
 
+    vLocalGear = in.Tb2l * vWhlBodyVec; // Get local frame wheel location
     gearLoc = in.Location.LocalToLocation(vLocalGear);
+
     // Compute the height of the theoretical location of the wheel (if strut is
     // not compressed) with respect to the ground level
     double height = gearLoc.GetContactPoint(t, contact, normal, terrainVel, dummy);
-    vGroundNormal = in.Tec2b * normal;
 
-    // The height returned above is the AGL and is expressed in the Z direction
+    if (height < 0.0) {
-    // of the ECEF coordinate frame. We now need to transform this height in
-    // actual compression of the strut (BOGEY) of in the normal direction to the
-    // ground (STRUCTURE)
-    double normalZ = (in.Tec2l*normal)(eZ);
-    double LGearProj = -(mTGear.Transposed() * vGroundNormal)(eZ);
-
-    switch (eContactType) {
-    case ctBOGEY:
-      compressLength = LGearProj > 0.0 ? height * normalZ / LGearProj : 0.0;
-      break;
-    case ctSTRUCTURE:
-      compressLength = height * normalZ / DotProduct(normal, normal);
-      break;
-    }
-
-    if (compressLength > 0.00) {
       WOW = true;
+      vGroundNormal = in.Tec2b * normal;
+
+      // The height returned by GetGroundCallback() is the AGL and is expressed
+      // in the Z direction of the local coordinate frame. We now need to transform
+      // this height in actual compression of the strut (BOGEY) or in the normal
+      // direction to the ground (STRUCTURE)
+      double normalZ = (in.Tec2l*normal)(eZ);
+      double LGearProj = -(mTGear.Transposed() * vGroundNormal)(eZ);
+      FGColumnVector3 vWhlDisplVec;
 
       // The following equations use the vector to the tire contact patch
       // including the strut compression.
-      FGColumnVector3 vWhlDisplVec;
-
       switch(eContactType) {
       case ctBOGEY:
+        compressLength = LGearProj > 0.0 ? height * normalZ / LGearProj : 0.0;
         vWhlDisplVec = mTGear * FGColumnVector3(0., 0., -compressLength);
         break;
       case ctSTRUCTURE:
+        compressLength = height * normalZ / DotProduct(normal, normal);
         vWhlDisplVec = compressLength * vGroundNormal;
         break;
       }
 
-      FGColumnVector3 vWhlContactVec = in.vWhlBodyVec[GearNumber] + vWhlDisplVec;
+      FGColumnVector3 vWhlContactVec = vWhlBodyVec + vWhlDisplVec;
       vActingXYZn = vXYZn + Tb2s * vWhlDisplVec;
       FGColumnVector3 vBodyWhlVel = in.PQR * vWhlContactVec;
       vBodyWhlVel += in.UVW - in.Tec2b * terrainVel;
@@ -334,15 +300,16 @@ const FGColumnVector3& FGLGear::GetBodyForces(void)
 
       InitializeReporting();
       ComputeSteeringAngle();
-      ComputeGroundCoordSys();
+      ComputeGroundFrame();
 
-      vLocalWhlVel = Transform().Transposed() * vBodyWhlVel;
+      vGroundWhlVel = mT.Transposed() * vBodyWhlVel;
 
       if (fdmex->GetTrimStatus())
         compressSpeed = 0.0; // Steady state is sought during trimming
       else {
-        compressSpeed = -vLocalWhlVel(eX);
+        compressSpeed = -vGroundWhlVel(eZ);
-        if (eContactType == ctBOGEY) compressSpeed /= LGearProj;
+        if (eContactType == ctBOGEY)
+          compressSpeed /= LGearProj;
       }
 
       ComputeVerticalStrutForce();
@@ -366,16 +333,24 @@ const FGColumnVector3& FGLGear::GetBodyForces(void)
       WheelSlip = 0.0;
       StrutForce = 0.0;
 
+      LMultiplier[ftRoll].value = 0.0;
+      LMultiplier[ftSide].value = 0.0;
+      LMultiplier[ftDynamic].value = 0.0;
+
       // Let wheel spin down slowly
       vWhlVelVec(eX) -= 13.0 * in.TotalDeltaT;
       if (vWhlVelVec(eX) < 0.0) vWhlVelVec(eX) = 0.0;
 
       // Return to neutral position between 1.0 and 0.8 gear pos.
-      SteerAngle *= max(GetGearUnitPos()-0.8, 0.0)/0.2;
+      SteerAngle *= max(gearPos-0.8, 0.0)/0.2;
 
       ResetReporting();
     }
   }
+  else if (gearPos < 0.01) { // Gear UP
+    WOW = false;
+    vWhlVelVec.InitMatrix();
+  }
 
   if (!fdmex->GetTrimStatus()) {
     ReportTakeoffOrLanding();
@@ -392,60 +367,28 @@ const FGColumnVector3& FGLGear::GetBodyForces(void)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Build a local "ground" coordinate system defined by
-//  eX : normal to the ground
+//  eX : projection of the rolling direction on the ground
-//  eY : projection of the rolling direction on the ground
+//  eY : projection of the sliping direction on the ground
-//  eZ : projection of the sliping direction on the ground
+//  eZ : normal to the ground
 
-void FGLGear::ComputeGroundCoordSys(void)
+void FGLGear::ComputeGroundFrame(void)
 {
-  // Euler angles are built up to create a local frame to describe the forces
+  FGColumnVector3 roll = mTGear * FGColumnVector3(cos(SteerAngle), sin(SteerAngle), 0.);
-  // applied to the gear by the ground. Here pitch, yaw and roll do not have
+  FGColumnVector3 side = vGroundNormal * roll;
-  // any physical meaning. It is just a convenient notation.
-  // First, "pitch" and "yaw" are determined in order to align eX with the
-  // ground normal.
-  if (vGroundNormal(eZ) < -1.0)
-    vOrient(ePitch) = 0.5*M_PI;
-  else if (1.0 < vGroundNormal(eZ))
-    vOrient(ePitch) = -0.5*M_PI;
-  else
-    vOrient(ePitch) = asin(-vGroundNormal(eZ));
 
-  if (fabs(vOrient(ePitch)) == 0.5*M_PI)
+  roll -= DotProduct(roll, vGroundNormal) * vGroundNormal;
-    vOrient(eYaw) = 0.;
+  roll.Normalize();
-  else
+  side.Normalize();
-    vOrient(eYaw) = atan2(vGroundNormal(eY), vGroundNormal(eX));
-  
-  vOrient(eRoll) = 0.;
-  UpdateCustomTransformMatrix();
 
-  if (eContactType == ctBOGEY) {
+  mT(eX,eX) = roll(eX);
-    // In the case of a bogey, the third angle "roll" is used to align the axis eY and eZ
+  mT(eY,eX) = roll(eY);
-    // to the rolling and sliping direction respectively.
+  mT(eZ,eX) = roll(eZ);
-    FGColumnVector3 updatedRollingAxis = Transform().Transposed() * mTGear
+  mT(eX,eY) = side(eX);
-                                       * FGColumnVector3(-sin(SteerAngle), cos(SteerAngle), 0.);
+  mT(eY,eY) = side(eY);
-
+  mT(eZ,eY) = side(eZ);
-    vOrient(eRoll) = atan2(updatedRollingAxis(eY), -updatedRollingAxis(eZ));
+  mT(eX,eZ) = vGroundNormal(eX);
-    UpdateCustomTransformMatrix();
+  mT(eY,eZ) = vGroundNormal(eY);
-  }
+  mT(eZ,eZ) = vGroundNormal(eZ);
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-void FGLGear::ComputeRetractionState(void)
-{
-  double gearPos = GetGearUnitPos();
-  if (gearPos < 0.01) {
-    GearUp   = true;
-    WOW      = false;
-    GearDown = false;
-    vWhlVelVec.InitMatrix();
-  } else if (gearPos > 0.99) {
-    GearDown = true;
-    GearUp   = false;
-  } else {
-    GearUp   = false;
-    GearDown = false;
-  }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -454,8 +397,8 @@ void FGLGear::ComputeRetractionState(void)
 void FGLGear::ComputeSlipAngle(void)
 {
 // Check that the speed is non-null otherwise use the current angle
-  if (vLocalWhlVel.Magnitude(eY,eZ) > 1E-3)
+  if (vGroundWhlVel.Magnitude(eX,eY) > 1E-3)
-    WheelSlip = -atan2(vLocalWhlVel(eZ), fabs(vLocalWhlVel(eY)))*radtodeg;
+    WheelSlip = -atan2(vGroundWhlVel(eY), fabs(vGroundWhlVel(eX)))*radtodeg;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -585,34 +528,10 @@ void FGLGear::CrashDetect(void)
 
 void FGLGear::ComputeBrakeForceCoefficient(void)
 {
-  switch (eBrakeGrp) {
+  BrakeFCoeff = rollingFCoeff;
-  case bgLeft:
+
-    BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgLeft]) +
+  if (eBrakeGrp != bgNone)
-                     staticFCoeff * in.BrakePos[bgLeft] );
+    BrakeFCoeff += in.BrakePos[eBrakeGrp] * (staticFCoeff - rollingFCoeff);
-    break;
-  case bgRight:
-    BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgRight]) +
-                     staticFCoeff * in.BrakePos[bgRight] );
-    break;
-  case bgCenter:
-    BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgCenter]) +
-                     staticFCoeff * in.BrakePos[bgCenter] );
-    break;
-  case bgNose:
-    BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgCenter]) +
-                     staticFCoeff * in.BrakePos[bgCenter] );
-    break;
-  case bgTail:
-    BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgCenter]) +
-                     staticFCoeff * in.BrakePos[bgCenter] );
-    break;
-  case bgNone:
-    BrakeFCoeff =  rollingFCoeff;
-    break;
-  default:
-    cerr << "Improper brake group membership detected for this gear." << endl;
-    break;
-  }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -651,8 +570,10 @@ void FGLGear::ComputeVerticalStrutForce(void)
 
   if (compressSpeed >= 0.0) {
 
-    if (eDampType == dtLinear)   dampForce = -compressSpeed * bDamp;
+    if (eDampType == dtLinear)
-    else         dampForce = -compressSpeed * compressSpeed * bDamp;
+      dampForce = -compressSpeed * bDamp;
+    else
+      dampForce = -compressSpeed * compressSpeed * bDamp;
 
   } else {
 
@@ -669,10 +590,10 @@ void FGLGear::ComputeVerticalStrutForce(void)
   switch (eContactType) {
   case ctBOGEY:
     // Project back the strut force in the local coordinate frame of the ground
-    vFn(eX) = StrutForce / (mTGear.Transposed()*vGroundNormal)(eZ);
+    vFn(eZ) = StrutForce / (mTGear.Transposed()*vGroundNormal)(eZ);
     break;
   case ctSTRUCTURE:
-    vFn(eX) = -StrutForce;
+    vFn(eZ) = -StrutForce;
     break;
   }
 
@@ -683,7 +604,7 @@ void FGLGear::ComputeVerticalStrutForce(void)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-double FGLGear::GetGearUnitPos(void)
+double FGLGear::GetGearUnitPos(void) const
 {
   // hack to provide backward compatibility to gear/gear-pos-norm property
   if( useFCSGearPos || in.FCSGearPos != 1.0 ) {
@@ -702,18 +623,19 @@ void FGLGear::ComputeJacobian(const FGColumnVector3& vWhlContactVec)
   // When the point of contact is moving, dynamic friction is used
   // This type of friction is limited to ctSTRUCTURE elements because their
   // friction coefficient is the same in every directions
-  if ((eContactType == ctSTRUCTURE) && (vLocalWhlVel.Magnitude(eY,eZ) > 1E-3)) {
+  if ((eContactType == ctSTRUCTURE) && (vGroundWhlVel.Magnitude(eX,eY) > 1E-3)) {
-    FGColumnVector3 velocityDirection = vLocalWhlVel;
+
+    FGColumnVector3 velocityDirection = vGroundWhlVel;
 
     StaticFriction = false;
 
-    velocityDirection(eX) = 0.;
+    velocityDirection(eZ) = 0.;
     velocityDirection.Normalize();
 
-    LMultiplier[ftDynamic].ForceJacobian = Transform()*velocityDirection;
+    LMultiplier[ftDynamic].ForceJacobian = mT * velocityDirection;
     LMultiplier[ftDynamic].MomentJacobian = vWhlContactVec * LMultiplier[ftDynamic].ForceJacobian;
     LMultiplier[ftDynamic].Max = 0.;
-    LMultiplier[ftDynamic].Min = -fabs(dynamicFCoeff * vFn(eX));
+    LMultiplier[ftDynamic].Min = -fabs(dynamicFCoeff * vFn(eZ));
 
     // The Lagrange multiplier value obtained from the previous iteration is kept
     // This is supposed to accelerate the convergence of the projected Gauss-Seidel
@@ -730,19 +652,19 @@ void FGLGear::ComputeJacobian(const FGColumnVector3& vWhlContactVec)
     // This cannot be handled properly by the so-called "dynamic friction".
     StaticFriction = true;
 
-    LMultiplier[ftRoll].ForceJacobian = Transform()*FGColumnVector3(0.,1.,0.);
+    LMultiplier[ftRoll].ForceJacobian = mT * FGColumnVector3(1.,0.,0.);
-    LMultiplier[ftSide].ForceJacobian = Transform()*FGColumnVector3(0.,0.,1.);
+    LMultiplier[ftSide].ForceJacobian = mT * FGColumnVector3(0.,1.,0.);
     LMultiplier[ftRoll].MomentJacobian = vWhlContactVec * LMultiplier[ftRoll].ForceJacobian;
     LMultiplier[ftSide].MomentJacobian = vWhlContactVec * LMultiplier[ftSide].ForceJacobian;
 
     switch(eContactType) {
     case ctBOGEY:
-      LMultiplier[ftRoll].Max = fabs(BrakeFCoeff * vFn(eX));
+      LMultiplier[ftRoll].Max = fabs(BrakeFCoeff * vFn(eZ));
-      LMultiplier[ftSide].Max = fabs(FCoeff * vFn(eX));
+      LMultiplier[ftSide].Max = fabs(FCoeff * vFn(eZ));
       break;
     case ctSTRUCTURE:
-      LMultiplier[ftRoll].Max = fabs(staticFCoeff * vFn(eX));
+      LMultiplier[ftRoll].Max = fabs(staticFCoeff * vFn(eZ));
-      LMultiplier[ftSide].Max = fabs(staticFCoeff * vFn(eX));
+      LMultiplier[ftSide].Max = LMultiplier[ftRoll].Max;
       break;
     }
 
@@ -768,11 +690,14 @@ void FGLGear::ComputeJacobian(const FGColumnVector3& vWhlContactVec)
 void FGLGear::UpdateForces(void)
 {
   if (StaticFriction) {
-    vFn(eY) = LMultiplier[ftRoll].value;
+    vFn(eX) = LMultiplier[ftRoll].value;
-    vFn(eZ) = LMultiplier[ftSide].value;
+    vFn(eY) = LMultiplier[ftSide].value;
+  }
+  else {
+    FGColumnVector3 forceDir = mT.Transposed() * LMultiplier[ftDynamic].ForceJacobian;
+    vFn(eX) = LMultiplier[ftDynamic].value * forceDir(eX);
+    vFn(eY) = LMultiplier[ftDynamic].value * forceDir(eY);
   }
-  else
-    vFn += LMultiplier[ftDynamic].value * (Transform ().Transposed() * LMultiplier[ftDynamic].ForceJacobian);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -892,11 +817,15 @@ void FGLGear::Report(ReportType repType)
 
 void FGLGear::Debug(int from)
 {
+  static const char* sSteerType[] = {"STEERABLE", "FIXED", "CASTERED" };
+  static const char* sBrakeGroup[] = {"NONE", "LEFT", "RIGHT", "CENTER", "NOSE", "TAIL"};
+  static const char* sContactType[] = {"BOGEY", "STRUCTURE" };
+
   if (debug_lvl <= 0) return;
 
   if (debug_lvl & 1) { // Standard console startup message output
     if (from == 0) { // Constructor - loading and initialization
-      cout << "    " << sContactType << " " << name          << endl;
+      cout << "    " << sContactType[eContactType] << " " << name          << endl;
       cout << "      Location: "         << vXYZn          << endl;
       cout << "      Spring Constant:  " << kSpring       << endl;
 
@@ -907,15 +836,15 @@ void FGLGear::Debug(int from)
 
       if (eDampTypeRebound == dtLinear)
         cout << "      Rebound Damping Constant: " << bDampRebound << " (linear)" << endl;
-      else 
+      else
         cout << "      Rebound Damping Constant: " << bDampRebound << " (square law)" << endl;
 
       cout << "      Dynamic Friction: " << dynamicFCoeff << endl;
       cout << "      Static Friction:  " << staticFCoeff  << endl;
       if (eContactType == ctBOGEY) {
         cout << "      Rolling Friction: " << rollingFCoeff << endl;
-        cout << "      Steering Type:    " << sSteerType    << endl;
+        cout << "      Steering Type:    " << sSteerType[eSteerType] << endl;
-        cout << "      Grouping:         " << sBrakeGroup   << endl;
+        cout << "      Grouping:         " << sBrakeGroup[eBrakeGrp] << endl;
         cout << "      Max Steer Angle:  " << maxSteerAngle << endl;
         cout << "      Retractable:      " << isRetractable  << endl;
       }
@@ -940,4 +869,3 @@ void FGLGear::Debug(int from)
 }
 
 } // namespace JSBSim
-

src/FDM/JSBSim/models/FGLGear.h

@@ -42,14 +42,13 @@ INCLUDES
 
 #include "models/propulsion/FGForce.h"
 #include "math/FGColumnVector3.h"
-#include "math/FGMatrix33.h"
 #include "math/LagrangeMultiplier.h"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_LGEAR "$Id: FGLGear.h,v 1.48 2011/10/31 14:54:41 bcoconni Exp $"
+#define ID_LGEAR "$Id: FGLGear.h,v 1.54 2012/04/01 17:05:51 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -178,7 +177,7 @@ CLASS DOCUMENTATION
         </contact>
 @endcode
     @author Jon S. Berndt
-    @version $Id: FGLGear.h,v 1.48 2011/10/31 14:54:41 bcoconni Exp $
+    @version $Id: FGLGear.h,v 1.54 2012/04/01 17:05:51 bcoconni Exp $
     @see Richard E. McFarland, "A Standard Kinematic Model for Flight Simulation at
      NASA-Ames", NASA CR-2497, January 1975
     @see Barnes W. McCormick, "Aerodynamics, Aeronautics, and Flight Mechanics",
@@ -208,10 +207,10 @@ public:
     FGMatrix33 Tec2b;
     FGColumnVector3 PQR;
     FGColumnVector3 UVW;
+    FGColumnVector3 vXYZcg; // CG coordinates expressed in the structural frame
     FGLocation Location;
     std::vector <double> SteerPosDeg;
     std::vector <double> BrakePos;
-    std::vector <FGColumnVector3> vWhlBodyVec;
     double FCSGearPos;
     double EmptyWeight;
   };
@@ -241,8 +240,13 @@ public:
   const FGColumnVector3& GetBodyForces(void);
 
   /// Gets the location of the gear in Body axes
-  const FGColumnVector3& GetBodyLocation(void) const { return in.vWhlBodyVec[GearNumber]; }
+  FGColumnVector3 GetBodyLocation(void) const {
-  double GetBodyLocation(int idx) const { return in.vWhlBodyVec[GearNumber](idx); }
+    return Ts2b * (vXYZn - in.vXYZcg);
+  }
+  double GetBodyLocation(int idx) const {
+    FGColumnVector3 vWhlBodyVec = Ts2b * (vXYZn - in.vXYZcg);
+    return vWhlBodyVec(idx);
+  }
 
   const FGColumnVector3& GetLocalGear(void) const { return vLocalGear; }
   double GetLocalGear(int idx) const { return vLocalGear(idx); }
@@ -257,15 +261,6 @@ public:
   double  GetCompVel(void) const {return compressSpeed; }
   /// Gets the gear compression force in pounds
   double  GetCompForce(void) const {return StrutForce;   }
-  double  GetBrakeFCoeff(void) const {return BrakeFCoeff;}
-
-  /// Gets the current normalized tire pressure
-  double  GetTirePressure(void) const { return TirePressureNorm; }
-  /// Sets the new normalized tire pressure
-  void    SetTirePressure(double p) { TirePressureNorm = p; }
-
-  /// Sets the brake value in percent (0 - 100)
-  void SetBrake(double bp) {brakePct = bp;}
 
   /// Sets the weight-on-wheels flag.
   void SetWOW(bool wow) {WOW = wow;}
@@ -285,8 +280,9 @@ public:
 
   bool GetSteerable(void) const    { return eSteerType != stFixed; }
   bool GetRetractable(void) const  { return isRetractable;   }
-  bool GetGearUnitUp(void) const   { return GearUp;          }
+  bool GetGearUnitUp(void) const   { return isRetractable ? (GetGearUnitPos() < 0.01) : false; }
-  bool GetGearUnitDown(void) const { return GearDown;        }
+  bool GetGearUnitDown(void) const { return isRetractable ? (GetGearUnitPos() > 0.99) : true; }
+
   double GetWheelRollForce(void) {
     UpdateForces();
     FGColumnVector3 vForce = mTGear.Transposed() * FGForce::GetBodyForces();
@@ -302,7 +298,7 @@ public:
   double GetWheelSlipAngle(void) const { return WheelSlip;       }
   double GetWheelVel(int axis) const   { return vWhlVelVec(axis);}
   bool IsBogey(void) const             { return (eContactType == ctBOGEY);}
-  double GetGearUnitPos(void);
+  double GetGearUnitPos(void) const;
   double GetSteerAngleDeg(void) const { return radtodeg*SteerAngle; }
 
   const struct Inputs& in;
@@ -311,13 +307,11 @@ public:
 
 private:
   int GearNumber;
-  static const FGMatrix33 Tb2s;
+  static const FGMatrix33 Tb2s, Ts2b;
   FGMatrix33 mTGear;
-  FGColumnVector3 vGearOrient;
   FGColumnVector3 vLocalGear;
-  FGColumnVector3 vWhlVelVec, vLocalWhlVel;     // Velocity of this wheel
+  FGColumnVector3 vWhlVelVec, vGroundWhlVel;     // Velocity of this wheel
-  FGColumnVector3 normal, vGroundNormal;
+  FGColumnVector3 vGroundNormal;
-  FGLocation contact, gearLoc;
   FGTable *ForceY_Table;
   double SteerAngle;
   double kSpring;
@@ -327,7 +321,6 @@ private:
   double compressSpeed;
   double staticFCoeff, dynamicFCoeff, rollingFCoeff;
   double Stiffness, Shape, Peak, Curvature; // Pacejka factors
-  double brakePct;
   double BrakeFCoeff;
   double maxCompLen;
   double SinkRate;
@@ -339,9 +332,7 @@ private:
   double MaximumStrutTravel;
   double FCoeff;
   double WheelSlip;
-  double TirePressureNorm;
   double GearPos;
-  bool   useFCSGearPos;
   bool WOW;
   bool lastWOW;
   bool FirstContact;
@@ -350,15 +341,9 @@ private:
   bool TakeoffReported;
   bool ReportEnable;
   bool isRetractable;
-  bool GearUp, GearDown;
-  bool Servicable;
   bool Castered;
   bool StaticFriction;
   std::string name;
-  std::string sSteerType;
-  std::string sBrakeGroup;
-  std::string sRetractable;
-  std::string sContactType;
 
   BrakeGroup  eBrakeGrp;
   ContactType eContactType;
@@ -372,13 +357,14 @@ private:
   FGGroundReactions* GroundReactions;
   FGPropertyManager* PropertyManager;
 
-  void ComputeRetractionState(void);
+  mutable bool useFCSGearPos;
+
   void ComputeBrakeForceCoefficient(void);
   void ComputeSteeringAngle(void);
   void ComputeSlipAngle(void);
   void ComputeSideForceCoefficient(void);
   void ComputeVerticalStrutForce(void);
-  void ComputeGroundCoordSys(void);
+  void ComputeGroundFrame(void);
   void ComputeJacobian(const FGColumnVector3& vWhlContactVec);
   void UpdateForces(void);
   void CrashDetect(void);

src/FDM/JSBSim/models/FGOutput.cpp

@@ -77,7 +77,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGOutput.cpp,v 1.63 2011/10/10 02:33:34 jentron Exp $";
+static const char *IdSrc = "$Id: FGOutput.cpp,v 1.67 2012/04/27 12:14:56 jberndt Exp $";
 static const char *IdHdr = ID_OUTPUT;
 
 // (stolen from FGFS native_fdm.cxx)
@@ -321,7 +321,7 @@ void FGOutput::DelimitedOutput(const string& fname)
     }
     if (SubSystems & ssMoments) {
       outstream << delimeter;
-      outstream << "L_{Aero} (ft-lbs)" + delimeter + "M_{Aero} ( ft-lbs)" + delimeter + "N_{Aero} (ft-lbs)" + delimeter;
+      outstream << "L_{Aero} (ft-lbs)" + delimeter + "M_{Aero} (ft-lbs)" + delimeter + "N_{Aero} (ft-lbs)" + delimeter;
       outstream << "L_{Prop} (ft-lbs)" + delimeter + "M_{Prop} (ft-lbs)" + delimeter + "N_{Prop} (ft-lbs)" + delimeter;
       outstream << "L_{Gear} (ft-lbs)" + delimeter + "M_{Gear} (ft-lbs)" + delimeter + "N_{Gear} (ft-lbs)" + delimeter;
       outstream << "L_{ext} (ft-lbs)" + delimeter + "M_{ext} (ft-lbs)" + delimeter + "N_{ext} (ft-lbs)" + delimeter;
@@ -981,15 +981,32 @@ void FGOutput::SocketStatusOutput(const string& out_str)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+bool FGOutput::Load(int subSystems, std::string protocol, std::string  type, std::string port, std::string name, double outRate, std::vector<FGPropertyManager *> & outputProperties)
+{
+  SetType(type);
+  SetRate(outRate);
+  SubSystems = subSystems;
+  OutputProperties = outputProperties;
+
+  if (((Type == otCSV) || (Type == otTab)) && (name != "cout") && (name !="COUT"))
+    name = FDMExec->GetRootDir() + name;
+
+  if (!port.empty() && (Type == otSocket || Type == otFlightGear)) {
+    SetProtocol(protocol);
+    socket = new FGfdmSocket(name, atoi(port.c_str()), Protocol);
+  } else {
+    BaseFilename = Filename = name;
+  }
+
+  Debug(2);
+  return true;
+}
+
 bool FGOutput::Load(Element* element)
 {
-  string parameter="";
+  int subSystems = 0;
-  string name="";
-  double OutRate = 0.0;
-  unsigned int port;
   Element *property_element;
-
+  std::vector<FGPropertyManager *> outputProperties;
-  string separator = "/";
 
   if (!DirectivesFile.empty()) { // A directives filename from the command line overrides
     output_file_name = DirectivesFile;      // one found in the config file.
@@ -1003,52 +1020,42 @@ bool FGOutput::Load(Element* element)
 
   if (!document) return false;
 
-  SetType(document->GetAttributeValue("type"));
+  string type = document->GetAttributeValue("type");
-
+  string name = document->GetAttributeValue("name");
-  name = document->GetAttributeValue("name");
+  string port = document->GetAttributeValue("port");
-  if (((Type == otCSV) || (Type == otTab)) && (name != "cout") && (name !="COUT"))
+  string protocol = document->GetAttributeValue("protocol");
-    name = FDMExec->GetRootDir() + name;
-
-  Port = document->GetAttributeValue("port");
-  if (!Port.empty() && (Type == otSocket || Type == otFlightGear)) {
-    port = atoi(Port.c_str());
-    SetProtocol(document->GetAttributeValue("protocol"));
-    socket = new FGfdmSocket(name, port, Protocol);
-  } else {
-    BaseFilename = Filename = name;
-  }
   if (!document->GetAttributeValue("rate").empty()) {
-    OutRate = document->GetAttributeValueAsNumber("rate");
+    rate = document->GetAttributeValueAsNumber("rate");
   } else {
-    OutRate = 1;
+    rate = 1;
   }
 
   if (document->FindElementValue("simulation") == string("ON"))
-    SubSystems += ssSimulation;
+    subSystems += ssSimulation;
   if (document->FindElementValue("aerosurfaces") == string("ON"))
-    SubSystems += ssAerosurfaces;
+    subSystems += ssAerosurfaces;
   if (document->FindElementValue("rates") == string("ON"))
-    SubSystems += ssRates;
+    subSystems += ssRates;
   if (document->FindElementValue("velocities") == string("ON"))
-    SubSystems += ssVelocities;
+    subSystems += ssVelocities;
   if (document->FindElementValue("forces") == string("ON"))
-    SubSystems += ssForces;
+    subSystems += ssForces;
   if (document->FindElementValue("moments") == string("ON"))
-    SubSystems += ssMoments;
+    subSystems += ssMoments;
   if (document->FindElementValue("atmosphere") == string("ON"))
-    SubSystems += ssAtmosphere;
+    subSystems += ssAtmosphere;
   if (document->FindElementValue("massprops") == string("ON"))
-    SubSystems += ssMassProps;
+    subSystems += ssMassProps;
   if (document->FindElementValue("position") == string("ON"))
-    SubSystems += ssPropagate;
+    subSystems += ssPropagate;
-  if (document->FindElementValue("coefficients") == string("ON"))
+  if (document->FindElementValue("coefficients") == string("ON") || document->FindElementValue("aerodynamics") == string("ON"))
-    SubSystems += ssAeroFunctions;
+    subSystems += ssAeroFunctions;
   if (document->FindElementValue("ground_reactions") == string("ON"))
-    SubSystems += ssGroundReactions;
+    subSystems += ssGroundReactions;
   if (document->FindElementValue("fcs") == string("ON"))
-    SubSystems += ssFCS;
+    subSystems += ssFCS;
   if (document->FindElementValue("propulsion") == string("ON"))
-    SubSystems += ssPropulsion;
+    subSystems += ssPropulsion;
   property_element = document->FindElement("property");
   while (property_element) {
     string property_str = property_element->GetDataLine();
@@ -1059,18 +1066,16 @@ bool FGOutput::Load(Element* element)
            << "  not be logged. You should check your configuration file."
            << reset << endl;
     } else {
-      OutputProperties.push_back(node);
+      outputProperties.push_back(node);
     }
     property_element = document->FindNextElement("property");
   }
 
-  SetRate(OutRate);
+  return Load(subSystems, protocol, type, port, name, rate, outputProperties);
-
-  Debug(2);
-
-  return true;
 }
 
+
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 void FGOutput::SetRate(double rtHz)

src/FDM/JSBSim/models/FGOutput.h

@@ -51,7 +51,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_OUTPUT "$Id: FGOutput.h,v 1.24 2011/11/10 12:06:14 jberndt Exp $"
+#define ID_OUTPUT "$Id: FGOutput.h,v 1.25 2012/02/07 23:15:37 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -124,7 +124,7 @@ CLASS DOCUMENTATION
     propulsion       ON|OFF
 </pre>
     NOTE that Time is always output with the data.
-    @version $Id: FGOutput.h,v 1.24 2011/11/10 12:06:14 jberndt Exp $
+    @version $Id: FGOutput.h,v 1.25 2012/02/07 23:15:37 bcoconni Exp $
  */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -167,6 +167,9 @@ public:
   bool Toggle(void) {enabled = !enabled; return enabled;}
 
   bool Load(Element* el);
+  bool Load(int subSystems, std::string protocol, std::string type, std::string port, 
+                            std::string name, double outRate,
+                            std::vector<FGPropertyManager *> & outputProperties);
   string GetOutputFileName(void) const {return Filename;}
 
   /// Subsystem types for specifying which will be output in the FDM data logging

src/FDM/JSBSim/models/FGPropagate.cpp

@@ -48,7 +48,16 @@ COMMENTS, REFERENCES,  and NOTES
     Wiley & Sons, 1979 ISBN 0-471-03032-5
 [5] Bernard Etkin, "Dynamics of Flight, Stability and Control", Wiley & Sons,
     1982 ISBN 0-471-08936-2
-[6] Erin Catto, "Iterative Dynamics with Temporal Coherence", February 22, 2005
+[6] S. Buss, "Accurate and Efficient Simulation of Rigid Body Rotations",
+    Technical Report, Department of Mathematics, University of California,
+    San Diego, 1999
+[7] Barker L.E., Bowles R.L. and Williams L.H., "Development and Application of
+    a Local Linearization Algorithm for the Integration of Quaternion Rate
+    Equations in Real-Time Flight Simulation Problems", NASA TN D-7347,
+    December 1973
+[8] Phillips W.F, Hailey C.E and Gebert G.A, "Review of Attitude Representations
+    Used for Aircraft Kinematics", Journal Of Aircraft Vol. 38, No. 4,
+    July-August 2001
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 INCLUDES
@@ -68,7 +77,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGPropagate.cpp,v 1.100 2011/11/06 18:14:51 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGPropagate.cpp,v 1.105 2012/03/26 21:26:11 bcoconni Exp $";
 static const char *IdHdr = ID_PROPAGATE;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -81,7 +90,7 @@ FGPropagate::FGPropagate(FGFDMExec* fdmex)
 {
   Debug(0);
   Name = "FGPropagate";
- 
+
   vInertialVelocity.InitMatrix();
 
   /// These define the indices use to select the various integrators.
@@ -200,7 +209,7 @@ This propagation is done using the current state values
 and current derivatives. Based on these values we compute an approximation to the
 state values for (now + dt).
 
-In the code below, variables named beginning with a small "v" refer to a 
+In the code below, variables named beginning with a small "v" refer to a
 a column vector, variables named beginning with a "T" refer to a transformation
 matrix. ECEF refers to Earth Centered Earth Fixed. ECI refers to Earth Centered
 Inertial.
@@ -214,12 +223,10 @@ bool FGPropagate::Run(bool Holding)
 
   double dt = in.DeltaT * rate;  // The 'stepsize'
 
-  RunPreFunctions();
-
   // Propagate rotational / translational velocity, angular /translational position, respectively.
 
-  Integrate(VState.vPQRi,             in.vPQRidot,          VState.dqPQRidot,          dt, integrator_rotational_rate);
   Integrate(VState.qAttitudeECI,      in.vQtrndot,          VState.dqQtrndot,          dt, integrator_rotational_position);
+  Integrate(VState.vPQRi,             in.vPQRidot,          VState.dqPQRidot,          dt, integrator_rotational_rate);
   Integrate(VState.vInertialPosition, VState.vInertialVelocity, VState.dqInertialVelocity, dt, integrator_translational_position);
   Integrate(VState.vInertialVelocity, in.vUVWidot,          VState.dqUVWidot,          dt, integrator_translational_rate);
 
@@ -240,10 +247,7 @@ bool FGPropagate::Run(bool Holding)
   //    vLocation vector.
   UpdateLocationMatrices();
 
-  // 5. Normalize the ECI Attitude quaternion
+  // 5. Update the "Orientation-based" transformation matrices from the updated
-  VState.qAttitudeECI.Normalize();
-
-  // 6. Update the "Orientation-based" transformation matrices from the updated 
   //    orientation quaternion and vLocation vector.
   UpdateBodyMatrices();
 
@@ -261,8 +265,6 @@ bool FGPropagate::Run(bool Holding)
   // Compute vehicle velocity wrt ECEF frame, expressed in Local horizontal frame.
   vVel = Tb2l * VState.vUVW;
 
-  RunPostFunctions();
-
   Debug(2);
   return false;
 }
@@ -271,7 +273,7 @@ bool FGPropagate::Run(bool Holding)
   // Transform the velocity vector of the body relative to the origin (Earth
   // center) to be expressed in the inertial frame, and add the vehicle velocity
   // contribution due to the rotation of the planet.
-  // Reference: See Stevens and Lewis, "Aircraft Control and Simulation", 
+  // Reference: See Stevens and Lewis, "Aircraft Control and Simulation",
   //            Second edition (2004), eqn 1.5-16c (page 50)
 
 void FGPropagate::CalculateInertialVelocity(void)
@@ -313,6 +315,8 @@ void FGPropagate::Integrate( FGColumnVector3& Integrand,
     break;
   case eNone: // do nothing, freeze translational rate
     break;
+  default:
+    break;
   }
 }
 
@@ -338,9 +342,84 @@ void FGPropagate::Integrate( FGQuaternion& Integrand,
     break;
   case eAdamsBashforth4: Integrand += (1/24.0)*dt*(55.0*ValDot[0] - 59.0*ValDot[1] + 37.0*ValDot[2] - 9.0*ValDot[3]);
     break;
+  case eBuss1:
+    {
+      // This is the first order method as described in Samuel R. Buss paper[6].
+      // The formula from Buss' paper is transposed below to quaternions and is
+      // actually the exact solution of the quaternion differential equation
+      // qdot = 1/2*w*q when w is constant.
+      Integrand = Integrand * QExp(0.5 * dt * VState.vPQRi);
+    }
+    return; // No need to normalize since the quaternion exponential is always normal
+  case eBuss2:
+    {
+      // This is the 'augmented second-order method' from S.R. Buss paper [6].
+      // Unlike Runge-Kutta or Adams-Bashforth, it is a one-pass second-order
+      // method (see reference [6]).
+      FGColumnVector3 wi = VState.vPQRi;
+      FGColumnVector3 wdoti = in.vPQRidot;
+      FGColumnVector3 omega = wi + 0.5*dt*wdoti + dt*dt/12.*wdoti*wi;
+      Integrand = Integrand * QExp(0.5 * dt * omega);
+    }
+    return; // No need to normalize since the quaternion exponential is always normal
+  case eLocalLinearization:
+    {
+      // This is the local linearization algorithm of Barker et al. (see ref. [7])
+      // It is also a one-pass second-order method. The code below is based on the
+      // more compact formulation issued from equation (107) of ref. [8]. The
+      // constants C1, C2, C3 and C4 have the same value than those in ref. [7] pp. 11
+      FGColumnVector3 wi = 0.5 * VState.vPQRi;
+      FGColumnVector3 wdoti = 0.5 * in.vPQRidot;
+      double omegak2 = DotProduct(VState.vPQRi, VState.vPQRi);
+      double omegak = omegak2 > 1E-6 ? sqrt(omegak2) : 1E-6;
+      double rhok = 0.5 * dt * omegak;
+      double C1 = cos(rhok);
+      double C2 = 2.0 * sin(rhok) / omegak;
+      double C3 = 4.0 * (1.0 - C1) / (omegak*omegak);
+      double C4 = 4.0 * (dt - C2) / (omegak*omegak);
+      FGColumnVector3 Omega = C2*wi + C3*wdoti + C4*wi*wdoti;
+      FGQuaternion q;
+
+      q(1) = C1 - C4*DotProduct(wi, wdoti);
+      q(2) = Omega(eP);
+      q(3) = Omega(eQ);
+      q(4) = Omega(eR);
+
+      Integrand = Integrand * q;
+
+      /* Cross check with ref. [7] pp.11-12 formulas and code pp. 20
+      double pk = VState.vPQRi(eP);
+      double qk = VState.vPQRi(eQ);
+      double rk = VState.vPQRi(eR);
+      double pdotk = in.vPQRidot(eP);
+      double qdotk = in.vPQRidot(eQ);
+      double rdotk = in.vPQRidot(eR);
+      double Ap = -0.25 * (pk*pdotk + qk*qdotk + rk*rdotk);
+      double Bp = 0.25 * (pk*qdotk - qk*pdotk);
+      double Cp = 0.25 * (pdotk*rk - pk*rdotk);
+      double Dp = 0.25 * (qk*rdotk - qdotk*rk);
+      double C2p = sin(rhok) / omegak;
+      double C3p = 2.0 * (1.0 - cos(rhok)) / (omegak*omegak);
+      double H = C1 + C4 * Ap;
+      double G = -C2p*rk - C3p*rdotk + C4*Bp;
+      double J = C2p*qk + C3p*qdotk - C4*Cp;
+      double K = C2p*pk + C3p*pdotk - C4*Dp;
+
+      cout << "q:       " << q << endl;
+
+      // Warning! In the paper of Barker et al. the quaternion components are not
+      // ordered the same way as in JSBSim (see equations (2) and (3) of ref. [7]
+      // as well as the comment just below equation (3))
+      cout << "FORTRAN: " << H << " , " << K << " , " << J << " , " << -G << endl;*/
+    }
+    break; // The quaternion q is not normal so the normalization needs to be done.
   case eNone: // do nothing, freeze rotational rate
     break;
+  default:
+    break;
   }
+
+  Integrand.Normalize();
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -443,7 +522,6 @@ void FGPropagate::SetVState(const VehicleState& vstate)
 {
   //ToDo: Shouldn't all of these be set from the vstate vector passed in?
   VState.vLocation = vstate.vLocation;
-  VState.vLocation.SetEarthPositionAngle(vstate.vLocation.GetEPA());
   Ti2ec = VState.vLocation.GetTi2ec(); // useless ?
   Tec2i = Ti2ec.Transposed();
   UpdateLocationMatrices();
@@ -566,7 +644,7 @@ void FGPropagate::bind(void)
   PropertyManager->Tie("attitude/roll-rad", this, (int)ePhi, (PMF)&FGPropagate::GetEuler);
   PropertyManager->Tie("attitude/pitch-rad", this, (int)eTht, (PMF)&FGPropagate::GetEuler);
   PropertyManager->Tie("attitude/heading-true-rad", this, (int)ePsi, (PMF)&FGPropagate::GetEuler);
-  
+
   PropertyManager->Tie("simulation/integrator/rate/rotational", (int*)&integrator_rotational_rate);
   PropertyManager->Tie("simulation/integrator/rate/translational", (int*)&integrator_translational_rate);
   PropertyManager->Tie("simulation/integrator/position/rotational", (int*)&integrator_rotational_position);
@@ -608,7 +686,7 @@ void FGPropagate::Debug(int from)
   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
   }
   if (debug_lvl & 8 && from == 2) { // Runtime state variables
-    cout << endl << fgblue << highint << left 
+    cout << endl << fgblue << highint << left
          << "  Propagation Report (English units: ft, degrees) at simulation time " << FDMExec->GetSimTime() << " seconds"
          << reset << endl;
     cout << endl;

src/FDM/JSBSim/models/FGPropagate.h

@@ -49,7 +49,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_PROPAGATE "$Id: FGPropagate.h,v 1.67 2011/11/09 22:07:17 bcoconni Exp $"
+#define ID_PROPAGATE "$Id: FGPropagate.h,v 1.69 2012/04/29 13:27:51 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -69,18 +69,18 @@ CLASS DOCUMENTATION
     state of the vehicle given the forces and moments that act on it. The
     integration accounts for a rotating Earth.
 
-    Integration of rotational and translation position and rate can be 
+    Integration of rotational and translation position and rate can be
     customized as needed or frozen by the selection of no integrator. The
-    selection of which integrator to use is done through the setting of 
+    selection of which integrator to use is done through the setting of
     the associated property. There are four properties which can be set:
-    
+
     @code
     simulation/integrator/rate/rotational
     simulation/integrator/rate/translational
     simulation/integrator/position/rotational
     simulation/integrator/position/translational
     @endcode
-    
+
     Each of the integrators listed above can be set to one of the following values:
 
     @code
@@ -93,7 +93,7 @@ CLASS DOCUMENTATION
     @endcode
 
     @author Jon S. Berndt, Mathias Froehlich, Bertrand Coconnier
-    @version $Id: FGPropagate.h,v 1.67 2011/11/09 22:07:17 bcoconni Exp $
+    @version $Id: FGPropagate.h,v 1.69 2012/04/29 13:27:51 bcoconni Exp $
   */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -156,12 +156,13 @@ public:
 
   /// Destructor
   ~FGPropagate();
-  
+
   /// These define the indices use to select the various integrators.
-  enum eIntegrateType {eNone = 0, eRectEuler, eTrapezoidal, eAdamsBashforth2, eAdamsBashforth3, eAdamsBashforth4};
+  enum eIntegrateType {eNone = 0, eRectEuler, eTrapezoidal, eAdamsBashforth2,
+                       eAdamsBashforth3, eAdamsBashforth4, eBuss1, eBuss2, eLocalLinearization};
 
   /** Initializes the FGPropagate class after instantiation and prior to first execution.
-      The base class FGModel::InitModel is called first, initializing pointers to the 
+      The base class FGModel::InitModel is called first, initializing pointers to the
       other FGModel objects (and others).  */
   bool InitModel(void);
 
@@ -169,7 +170,7 @@ public:
 
   /** Runs the state propagation model; called by the Executive
       Can pass in a value indicating if the executive is directing the simulation to Hold.
-      @param Holding if true, the executive has been directed to hold the sim from 
+      @param Holding if true, the executive has been directed to hold the sim from
                      advancing time. Some models may ignore this flag, such as the Input
                      model, which may need to be active to listen on a socket for the
                      "Resume" command to be given.
@@ -188,7 +189,7 @@ public:
               expressed in Local horizontal frame.
   */
   const FGColumnVector3& GetVel(void) const { return vVel; }
-  
+
   /** Retrieves the body frame vehicle velocity vector.
       The vector returned is represented by an FGColumnVector reference. The vector
       for the velocity in Body frame is organized (Vx, Vy, Vz). The vector
@@ -200,7 +201,7 @@ public:
       @return The body frame vehicle velocity vector in ft/sec.
   */
   const FGColumnVector3& GetUVW(void) const { return VState.vUVW; }
-  
+
   /** Retrieves the body angular rates vector, relative to the ECEF frame.
       Retrieves the body angular rates (p, q, r), which are calculated by integration
       of the angular acceleration.
@@ -214,7 +215,7 @@ public:
       @return The body frame angular rates in rad/sec.
   */
   const FGColumnVector3& GetPQR(void) const {return VState.vPQR;}
-  
+
   /** Retrieves the body angular rates vector, relative to the ECI (inertial) frame.
       Retrieves the body angular rates (p, q, r), which are calculated by integration
       of the angular acceleration.
@@ -556,7 +557,6 @@ private:
   FGColumnVector3 vVel;
   FGColumnVector3 vInertialVelocity;
   FGColumnVector3 vLocation;
-  FGColumnVector3 vDeltaXYZEC;
   FGMatrix33 Tec2b;
   FGMatrix33 Tb2ec;
   FGMatrix33 Tl2b;   // local to body frame matrix copy for immediate local use
@@ -569,7 +569,7 @@ private:
   FGMatrix33 Tb2i;   // body to ECI frame rotation matrix
   FGMatrix33 Ti2l;
   FGMatrix33 Tl2i;
-  
+
   double VehicleRadius;
   FGColumnVector3 LocalTerrainVelocity, LocalTerrainAngularVelocity;
 

src/FDM/JSBSim/models/FGPropulsion.cpp

@@ -66,7 +66,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGPropulsion.cpp,v 1.52 2011/10/31 14:54:41 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGPropulsion.cpp,v 1.61 2012/04/14 18:10:44 bcoconni Exp $";
 static const char *IdHdr = ID_PROPULSION;
 
 extern short debug_lvl;
@@ -204,7 +204,7 @@ void FGPropulsion::ConsumeFuel(FGEngine* engine)
   unsigned int TanksWithOxidizer=0, CurrentOxidizerTankPriority=1;
   vector <int> FeedListFuel, FeedListOxi;
   bool Starved = true; // Initially set Starved to true. Set to false in code below.
-//  bool hasOxTanks = false;
+  bool hasOxTanks = false;
 
   // For this engine,
   // 1) Count how many fuel tanks with the current priority level have fuel
@@ -237,6 +237,9 @@ void FGPropulsion::ConsumeFuel(FGEngine* engine)
     if (TanksWithFuel == 0) CurrentFuelTankPriority++; // No tanks at this priority, try next priority
   }
 
+  bool FuelStarved = Starved;
+  Starved = true;
+
   // Process Oxidizer tanks, if any
   if (engine->GetType() == FGEngine::etRocket) {
     while ((TanksWithOxidizer == 0) && (CurrentOxidizerTankPriority <= numTanks)) {
@@ -250,7 +253,7 @@ void FGPropulsion::ConsumeFuel(FGEngine* engine)
             // Skip this here (done above)
             break;
           case FGTank::ttOXIDIZER:
-//            hasOxTanks = true;
+            hasOxTanks = true;
             if (Tank->GetContents() > 0.0 && Tank->GetSelected() && TankPriority == CurrentOxidizerTankPriority) {
               TanksWithOxidizer++;
               if (TanksWithFuel > 0) Starved = false;
@@ -264,10 +267,13 @@ void FGPropulsion::ConsumeFuel(FGEngine* engine)
     }
   }
 
-  engine->SetStarved(Starved); // Tanks can be refilled, so be sure to reset engine Starved flag here.
+  bool OxiStarved = Starved;
+
+  engine->SetStarved(FuelStarved || (hasOxTanks && OxiStarved)); // Tanks can be refilled, so be sure to reset engine Starved flag here.
 
   // No fuel or fuel/oxidizer found at any priority!
-  if (Starved) return;
+//  if (Starved) return;
+  if (FuelStarved || (hasOxTanks && OxiStarved)) return;
 
   double FuelToBurn = engine->CalcFuelNeed();            // How much fuel does this engine need?
   double FuelNeededPerTank = FuelToBurn / TanksWithFuel; // Determine fuel needed per tank.  
@@ -467,18 +473,18 @@ string FGPropulsion::FindEngineFullPathname(const string& engine_filename)
   fullpath = enginePath + separator;
   localpath = aircraftPath + separator + "Engines" + separator;
 
-  engine_file.open(string(fullpath + engine_filename + ".xml").c_str());
+  engine_file.open(string(localpath + engine_filename + ".xml").c_str());
   if ( !engine_file.is_open()) {
-    engine_file.open(string(localpath + engine_filename + ".xml").c_str());
+    engine_file.open(string(fullpath + engine_filename + ".xml").c_str());
       if ( !engine_file.is_open()) {
         cerr << " Could not open engine file: " << engine_filename << " in path "
              << fullpath << " or " << localpath << endl;
         return string("");
       } else {
-        return string(localpath + engine_filename + ".xml");
+        return string(fullpath + engine_filename + ".xml");
       }
   }
-  return string(fullpath + engine_filename + ".xml");
+  return string(localpath + engine_filename + ".xml");
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -495,12 +501,12 @@ ifstream* FGPropulsion::FindEngineFile(const string& engine_filename)
   fullpath = enginePath + separator;
   localpath = aircraftPath + separator + "Engines" + separator;
 
-  engine_file->open(string(fullpath + engine_filename + ".xml").c_str());
+  engine_file->open(string(localpath + engine_filename + ".xml").c_str());
   if ( !engine_file->is_open()) {
-    engine_file->open(string(localpath + engine_filename + ".xml").c_str());
+    engine_file->open(string(fullpath + engine_filename + ".xml").c_str());
       if ( !engine_file->is_open()) {
         cerr << " Could not open engine file: " << engine_filename << " in path "
-             << fullpath << " or " << localpath << endl;
+             << localpath << " or " << fullpath << endl;
       }
   }
   return engine_file;
@@ -527,6 +533,9 @@ string FGPropulsion::GetPropulsionStrings(const string& delimiter) const
     else if (Tanks[i]->GetType() == FGTank::ttOXIDIZER) buf << delimiter << "Oxidizer Tank " << i;
   }
 
+  PropulsionStrings += buf.str();
+  buf.str("");
+
   return PropulsionStrings;
 }
 
@@ -551,6 +560,9 @@ string FGPropulsion::GetPropulsionValues(const string& delimiter) const
     buf << Tanks[i]->GetContents();
   }
 
+  PropulsionValues += buf.str();
+  buf.str("");
+
   return PropulsionValues;
 }
 
@@ -588,9 +600,7 @@ const FGColumnVector3& FGPropulsion::GetTanksMoment(void)
 {
   vXYZtank_arm.InitMatrix();
   for (unsigned int i=0; i<Tanks.size(); i++) {
-    vXYZtank_arm(eX) += Tanks[i]->GetXYZ(eX) * Tanks[i]->GetContents();
+    vXYZtank_arm += Tanks[i]->GetXYZ() * Tanks[i]->GetContents();
-    vXYZtank_arm(eY) += Tanks[i]->GetXYZ(eY) * Tanks[i]->GetContents();
-    vXYZtank_arm(eZ) += Tanks[i]->GetXYZ(eZ) * Tanks[i]->GetContents();
   }
   return vXYZtank_arm;
 }
@@ -610,6 +620,7 @@ double FGPropulsion::GetTanksWeight(void) const
 
 const FGMatrix33& FGPropulsion::CalculateTankInertias(void)
 {
+  const FGMatrix33 Ts2b(-inchtoft, 0., 0., 0., inchtoft, 0., 0., 0., -inchtoft);
   unsigned int size;
 
   size = Tanks.size();
@@ -618,8 +629,10 @@ const FGMatrix33& FGPropulsion::CalculateTankInertias(void)
   tankJ = FGMatrix33();
 
   for (unsigned int i=0; i<size; i++) {
+    FGColumnVector3 vTankBodyVec = Ts2b * (in.vXYZcg - Tanks[i]->GetXYZ());
+
     tankJ += FDMExec->GetMassBalance()->GetPointmassInertia( lbtoslug * Tanks[i]->GetContents(),
-                                               Tanks[i]->GetXYZ() );
+                                                             vTankBodyVec);
     tankJ(1,1) += Tanks[i]->GetIxx();
     tankJ(2,2) += Tanks[i]->GetIyy();
     tankJ(3,3) += Tanks[i]->GetIzz();

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

@@ -2,7 +2,7 @@
 
  Module:       FGMSIS.cpp
  Author:       David Culp
-               (incorporated into C++ JSBSim class heirarchy, see model authors below)
+               (incorporated into C++ JSBSim class hierarchy, see model authors below)
  Date started: 12/14/03
  Purpose:      Models the MSIS-00 atmosphere
 
@@ -66,7 +66,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGMSIS.cpp,v 1.18 2011/06/21 13:54:40 jberndt Exp $";
+static const char *IdSrc = "$Id: FGMSIS.cpp,v 1.19 2011/12/11 17:03:05 bcoconni Exp $";
 static const char *IdHdr = ID_MSIS;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -161,12 +161,10 @@ bool MSIS::Run(bool Holding)
   if (FGModel::Run(Holding)) return true;
   if (Holding) return false;
 
-  RunPreFunctions();
-
   double h = FDMExec->GetPropagate()->GetAltitudeASL();
 
   //do temp, pressure, and density first
-//  if (!useExternal) {
+  //if (!useExternal) {
     // get sea-level values
     Calculate(FDMExec->GetAuxiliary()->GetDayOfYear(),
               FDMExec->GetAuxiliary()->GetSecondsInDay(),
@@ -188,14 +186,12 @@ bool MSIS::Run(bool Holding)
               h,
               FDMExec->GetPropagate()->GetLocation().GetLatitudeDeg(),
               FDMExec->GetPropagate()->GetLocation().GetLongitudeDeg());
-//    intTemperature = output.t[1] * 1.8;
+    //intTemperature = output.t[1] * 1.8;
-//    intDensity     = output.d[5] * 1.940321;
+    //intDensity     = output.d[5] * 1.940321;
-//    intPressure    = 1716.488 * intDensity * intTemperature;
+    //intPressure    = 1716.488 * intDensity * intTemperature;
     //cout << "T=" << intTemperature << " D=" << intDensity << " P=";
     //cout << intPressure << " a=" << soundspeed << endl;
-//  }
+  //}
-
-  RunPostFunctions();
 
   Debug(2);
 

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

@@ -50,7 +50,7 @@ INCLUDES
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGStandardAtmosphere.cpp,v 1.20 2011/09/18 12:06:21 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGStandardAtmosphere.cpp,v 1.21 2012/04/13 13:18:27 jberndt Exp $";
 static const char *IdHdr = ID_STANDARDATMOSPHERE;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -170,7 +170,7 @@ double FGStandardAtmosphere::GetPressure(double altitude) const
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-void FGStandardAtmosphere::SetPressureSL(double pressure, ePressure unit)
+void FGStandardAtmosphere::SetPressureSL(ePressure unit, double pressure)
 {
   double press = ConvertToPSF(pressure, unit);
 
@@ -421,6 +421,9 @@ void FGStandardAtmosphere::bind(void)
   PropertyManager->Tie("atmosphere/SL-graded-delta-T", this, eRankine,
                                     (PMFi)&FGStandardAtmosphere::GetTemperatureDeltaGradient,
                                     (PMF)&FGStandardAtmosphere::SetSLTemperatureGradedDelta);
+  PropertyManager->Tie("atmosphere/P-sl-psf", this, ePSF,
+                                   (PMFi)&FGStandardAtmosphere::GetPressureSL,
+                                   (PMF)&FGStandardAtmosphere::SetPressureSL);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

@@ -46,7 +46,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_STANDARDATMOSPHERE "$Id: FGStandardAtmosphere.h,v 1.16 2011/09/18 12:06:21 bcoconni Exp $"
+#define ID_STANDARDATMOSPHERE "$Id: FGStandardAtmosphere.h,v 1.17 2012/04/13 13:18:27 jberndt Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -93,7 +93,7 @@ consistently and accurately calculated.
 
   @author Jon Berndt
   @see "U.S. Standard Atmosphere, 1976", NASA TM-X-74335
-  @version $Id: FGStandardAtmosphere.h,v 1.16 2011/09/18 12:06:21 bcoconni Exp $
+  @version $Id: FGStandardAtmosphere.h,v 1.17 2012/04/13 13:18:27 jberndt Exp $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -226,10 +226,10 @@ public:
   /** Sets the sea level pressure for modeling an off-standard pressure
       profile. This could be useful in the case where the pressure at an
       airfield is known or set for a particular simulation run.
-      @param pressure The pressure in the units specified (PSF by default).
+      @param pressure The pressure in the units specified.
       @param unit the unit of measure that the specified pressure is
                        supplied in.*/
-  virtual void SetPressureSL(double pressure, ePressure unit=ePSF);
+  virtual void SetPressureSL(ePressure unit, double pressure);
 
   /** Resets the sea level to the Standard sea level pressure, and recalculates
       dependent parameters so that the pressure calculations are standard. */

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

@@ -51,7 +51,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGWinds.cpp,v 1.6 2011/11/10 12:02:34 jberndt Exp $";
+static const char *IdSrc = "$Id: FGWinds.cpp,v 1.7 2011/12/11 17:03:05 bcoconni Exp $";
 static const char *IdHdr = ID_WINDS;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -130,8 +130,6 @@ bool FGWinds::Run(bool Holding)
   if (FGModel::Run(Holding)) return true;
   if (Holding) return false;
 
-  RunPreFunctions();
-
   if (turbType != ttNone) Turbulence(in.AltitudeASL);
   if (oneMinusCosineGust.gustProfile.Running) CosineGust();
 
@@ -141,8 +139,6 @@ bool FGWinds::Run(bool Holding)
   if (vWindNED(eX) != 0.0) psiw = atan2( vWindNED(eY), vWindNED(eX) );
   if (psiw < 0) psiw += 2*M_PI;
 
-  RunPostFunctions();
-
   Debug(2);
   return false;
 }
@@ -178,7 +174,7 @@ void FGWinds::SetWindPsi(double dir)
 {
   double mag = GetWindspeed();
   psiw = dir;
-  SetWindspeed(mag);  
+  SetWindspeed(mag);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -187,11 +183,11 @@ void FGWinds::Turbulence(double h)
 {
   switch (turbType) {
 
-  case ttCulp: { 
+  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;
@@ -212,7 +208,7 @@ void FGWinds::Turbulence(double h)
     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;
@@ -225,7 +221,7 @@ void FGWinds::Turbulence(double h)
     vTurbulenceNED(3)+= delta;
     if (in.DistanceAGL/in.wingspan < 3.0)
         vTurbulenceNED(3) *= in.DistanceAGL/in.wingspan * 0.3333;
- 
+
     // Yaw component of turbulence.
     vTurbulenceNED(1) = sin( delta * 3.0 );
     vTurbulenceNED(2) = cos( delta * 3.0 );
@@ -501,7 +497,7 @@ void FGWinds::bind(void)
                                                           (PMFd)&FGWinds::SetGustNED);
   PropertyManager->Tie("atmosphere/gust-down-fps",  this, eDown, (PMF)&FGWinds::GetGustNED,
                                                           (PMFd)&FGWinds::SetGustNED);
-  
+
   // User-specified 1 - cosine gust parameters (in specified frame)
   PropertyManager->Tie("atmosphere/cosine-gust/startup-duration-sec", this, (Ptr)0L, &FGWinds::StartupGustDuration);
   PropertyManager->Tie("atmosphere/cosine-gust/steady-duration-sec", this, (Ptr)0L, &FGWinds::SteadyGustDuration);
@@ -592,7 +588,7 @@ void FGWinds::Debug(int from)
   }
   if (debug_lvl & 16) { // Sanity checking
   }
-  if (debug_lvl & 128) { // 
+  if (debug_lvl & 128) { //
   }
   if (debug_lvl & 64) {
     if (from == 0) { // Constructor

src/FDM/JSBSim/models/flight_control/FGAccelerometer.h

@@ -47,7 +47,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_ACCELEROMETER "$Id: FGAccelerometer.h,v 1.5 2011/07/17 13:51:23 jberndt Exp $"
+#define ID_ACCELEROMETER "$Id: FGAccelerometer.h,v 1.6 2012/01/08 12:39:14 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -71,7 +71,17 @@ Syntax:
 
 @code
 <accelerometer name="name">
-  <input> property </input>
+  <location unit="{IN | M}">
+    <x> number </x>
+    <y> number </y>
+    <z> number </z>
+  </location>
+  <orientation unit="{RAD | DEG}">
+    <pitch> {number} </pitch>
+    <roll> {number} </roll>
+    <yaw> {number} </yaw>
+  </orientation>
+  <axis> {X | Y | Z} </axis>
   <lag> number </lag>
   <noise variation="PERCENT|ABSOLUTE"> number </noise>
   <quantization name="name">
@@ -80,6 +90,7 @@ Syntax:
     <max> number </max>
   </quantization>
   <drift_rate> number </drift_rate>
+  <gain> number </gain>
   <bias> number </bias>
 </accelerometer>
 @endcode
@@ -87,11 +98,16 @@ Syntax:
 Example:
 
 @code
-<accelerometer name="aero/accelerometer/qbar">
+<accelerometer name="aero/accelerometer/right_tip_wing">
-  <input> aero/qbar </input>
+  <location unit="IN">
+    <x> 43.2 </x>
+    <y> 214. </y>
+    <z> 59.4 </z>
+  </location>
+  <axis> Z </axis>
   <lag> 0.5 </lag>
   <noise variation="PERCENT"> 2 </noise>
-  <quantization name="aero/accelerometer/quantized/qbar">
+  <quantization name="aero/accelerometer/quantized/right_tip_wing">
     <bits> 12 </bits>
     <min> 0 </min>
     <max> 400 </max>
@@ -111,7 +127,7 @@ even varying all the way from 0.95 to 1.05 in adjacent frames - whatever the del
 time.
 
 @author Jon S. Berndt
-@version $Revision: 1.5 $
+@version $Revision: 1.6 $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/FDM/JSBSim/models/flight_control/FGActuator.cpp

@@ -43,7 +43,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGActuator.cpp,v 1.22 2011/07/12 21:40:32 jentron Exp $";
+static const char *IdSrc = "$Id: FGActuator.cpp,v 1.23 2012/04/08 15:04:41 jberndt Exp $";
 static const char *IdHdr = ID_ACTUATOR;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -66,6 +66,7 @@ FGActuator::FGActuator(FGFCS* fcs, Element* element) : FGFCSComponent(fcs, eleme
   fail_zero = fail_hardover = fail_stuck = false;
   ca = cb = 0.0;
   initialized = 0;
+  saturated = false;
 
   if ( element->FindElement("deadband_width") ) {
     deadband_width = element->FindElementValueAsNumber("deadband_width");
@@ -132,6 +133,13 @@ bool FGActuator::Run(void )
   initialized = 1;
 
   Clip();
+
+  if (clip) {
+    saturated = false;
+    if (Output >= clipmax) saturated = true;
+    else if (Output <= clipmin) saturated = true;
+  }
+
   if (IsOutput) SetOutput();
 
   return true;
@@ -225,10 +233,12 @@ void FGActuator::bind(void)
   const string tmp_zero = tmp + "/malfunction/fail_zero";
   const string tmp_hardover = tmp + "/malfunction/fail_hardover";
   const string tmp_stuck = tmp + "/malfunction/fail_stuck";
+  const string tmp_sat = tmp + "/saturated";
 
   PropertyManager->Tie( tmp_zero, this, &FGActuator::GetFailZero, &FGActuator::SetFailZero);
   PropertyManager->Tie( tmp_hardover, this, &FGActuator::GetFailHardover, &FGActuator::SetFailHardover);
   PropertyManager->Tie( tmp_stuck, this, &FGActuator::GetFailStuck, &FGActuator::SetFailStuck);
+  PropertyManager->Tie( tmp_sat, this, &FGActuator::IsSaturated);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/FDM/JSBSim/models/flight_control/FGActuator.h

@@ -44,7 +44,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_ACTUATOR "$Id: FGActuator.h,v 1.12 2011/07/12 21:40:32 jentron Exp $"
+#define ID_ACTUATOR "$Id: FGActuator.h,v 1.13 2012/04/08 15:04:41 jberndt Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -111,7 +111,7 @@ Example:
 @endcode
 
 @author Jon S. Berndt
-@version $Revision: 1.12 $
+@version $Revision: 1.13 $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -135,13 +135,14 @@ public:
   /** This function fails the actuator to zero. The motion to zero
       will flow through the lag, hysteresis, and rate limiting
       functions if those are activated. */
-  inline void SetFailZero(bool set) {fail_zero = set;}
+  void SetFailZero(bool set) {fail_zero = set;}
-  inline void SetFailHardover(bool set) {fail_hardover = set;}
+  void SetFailHardover(bool set) {fail_hardover = set;}
-  inline void SetFailStuck(bool set) {fail_stuck = set;}
+  void SetFailStuck(bool set) {fail_stuck = set;}
 
-  inline bool GetFailZero(void) const {return fail_zero;}
+  bool GetFailZero(void) const {return fail_zero;}
-  inline bool GetFailHardover(void) const {return fail_hardover;}
+  bool GetFailHardover(void) const {return fail_hardover;}
-  inline bool GetFailStuck(void) const {return fail_stuck;}
+  bool GetFailStuck(void) const {return fail_stuck;}
+  bool IsSaturated(void) const {return saturated;}
   
 private:
   double span;
@@ -161,6 +162,7 @@ private:
   bool fail_hardover;
   bool fail_stuck;
   bool initialized;
+  bool saturated;
 
   void Hysteresis(void);
   void Lag(void);

src/FDM/JSBSim/models/flight_control/FGPID.cpp

@@ -44,7 +44,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGPID.cpp,v 1.19 2011/05/05 11:44:11 jberndt Exp $";
+static const char *IdSrc = "$Id: FGPID.cpp,v 1.20 2012/05/10 12:10:48 jberndt Exp $";
 static const char *IdHdr = ID_PID;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -65,6 +65,13 @@ FGPID::FGPID(FGFCS* fcs, Element* element) : FGFCSComponent(fcs, element)
   I_out_total = 0.0;
   Input_prev = Input_prev2 = 0.0;
   Trigger = 0;
+  ProcessVariableDot = 0;
+  IsStandard = false;
+  IntType = eNone;       // No integrator initially defined.
+
+  string pid_type = element->GetAttributeValue("type");
+
+  if (pid_type == "standard") IsStandard = true;
 
   if ( element->FindElement("kp") ) {
     kp_string = element->FindElementValue("kp");
@@ -81,6 +88,20 @@ FGPID::FGPID(FGFCS* fcs, Element* element) : FGFCSComponent(fcs, element)
 
   if ( element->FindElement("ki") ) {
     ki_string = element->FindElementValue("ki");
+
+    string integ_type = element->FindElement("ki")->GetAttributeValue("type");
+    if (integ_type == "rect") {            // Use rectangular integration
+      IntType = eRectEuler;
+    } else if (integ_type == "trap") {     // Use trapezoidal integration
+      IntType = eTrapezoidal;
+    } else if (integ_type == "ab2") {      // Use Adams Bashforth 2nd order integration
+      IntType = eAdamsBashforth2;
+    } else if (integ_type == "ab3") {      // Use Adams Bashforth 3rd order integration
+      IntType = eAdamsBashforth3;
+    } else {                               // Use default Adams Bashforth 2nd order integration
+      IntType = eAdamsBashforth2;
+    }
+
     if (!is_number(ki_string)) { // property
       if (ki_string[0] == '-') {
        KiPropertySign = -1.0;
@@ -105,6 +126,10 @@ FGPID::FGPID(FGFCS* fcs, Element* element) : FGFCSComponent(fcs, element)
     }
   }
 
+  if (element->FindElement("pvdot")) {
+    ProcessVariableDot =  PropertyManager->GetNode(element->FindElementValue("pvdot"));
+  }
+
   if (element->FindElement("trigger")) {
     Trigger =  PropertyManager->GetNode(element->FindElementValue("trigger"));
   }
@@ -126,7 +151,7 @@ FGPID::~FGPID()
 bool FGPID::Run(void )
 {
   double I_out_delta = 0.0;
-  double P_out, D_out;
+  double Dval = 0;
 
   Input = InputNodes[0]->getDoubleValue() * InputSigns[0];
 
@@ -134,30 +159,51 @@ bool FGPID::Run(void )
   if (KiPropertyNode != 0) Ki = KiPropertyNode->getDoubleValue() * KiPropertySign;
   if (KdPropertyNode != 0) Kd = KdPropertyNode->getDoubleValue() * KdPropertySign;
 
-  P_out = Kp * Input;
+  if (ProcessVariableDot) {
-  D_out = (Kd / dt) * (Input - Input_prev);
+    Dval = ProcessVariableDot->getDoubleValue();
+  } else {
+    Dval = (Input - Input_prev)/dt;
+  }
 
   // Do not continue to integrate the input to the integrator if a wind-up
   // condition is sensed - that is, if the property pointed to by the trigger
   // element is non-zero. Reset the integrator to 0.0 if the Trigger value
   // is negative.
 
-  if (Trigger != 0) {
+  double test = 0.0;
-    double test = Trigger->getDoubleValue();
+  if (Trigger != 0) test = Trigger->getDoubleValue();
-    if (fabs(test) < 0.000001) {
+
-      // I_out_delta = Ki * dt * Input;                         // Normal rectangular integrator
+  if (fabs(test) < 0.000001) {
+    switch(IntType) {
+    case eRectEuler:
+      I_out_delta = Ki * dt * Input;                         // Normal rectangular integrator
+      break;
+    case eTrapezoidal:
+      I_out_delta = (Ki/2.0) * dt * (Input + Input_prev);    // Trapezoidal integrator
+      break;
+    case eAdamsBashforth2:
       I_out_delta = Ki * dt * (1.5*Input - 0.5*Input_prev);  // 2nd order Adams Bashforth integrator
+      break;
+    case eAdamsBashforth3:                                   // 3rd order Adams Bashforth integrator
+      I_out_delta = (Ki/12.0) * dt * (23.0*Input - 16.0*Input_prev + 5.0*Input_prev2);
+      break;
+    case eNone:
+      // No integator is defined or used.
+      I_out_delta = 0.0;
+      break;
     }
-    if (test < 0.0)            I_out_total = 0.0;  // Reset integrator to 0.0
-  } else { // no anti-wind-up trigger defined
-    // I_out_delta = Ki * dt * Input;
-    I_out_delta = Ki * dt * (1.5*Input - 0.5*Input_prev);  // 2nd order Adams Bashforth integrator
   }
+
+  if (test < 0.0) I_out_total = 0.0;  // Reset integrator to 0.0
   
   I_out_total += I_out_delta;
 
-  Output = P_out + I_out_total + D_out;
+  if (IsStandard) {
-  
+    Output = Kp * (Input + I_out_total + Kd*Dval);
+  } else {
+    Output = Kp*Input + I_out_total + Kd*Dval;
+  }
+
   Input_prev = Input;
   Input_prev2 = Input_prev;
 

src/FDM/JSBSim/models/flight_control/FGPID.h

@@ -44,7 +44,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_PID "$Id: FGPID.h,v 1.12 2009/10/24 22:59:30 jberndt Exp $"
+#define ID_PID "$Id: FGPID.h,v 1.13 2012/05/10 12:10:48 jberndt Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -64,26 +64,59 @@ CLASS DOCUMENTATION
 <h3>Configuration Format:</h3>
 
 @code
-<pid name="{string}">
+<pid name="{string}" [type="standard"]>
-  <kp> {number} </kp>
+  <kp> {number|property} </kp>
-  <ki> {number} </ki>
+  <ki> {number|property} </ki>
-  <kd> {number} </kd>
+  <kd> {number|property} </kd>
-  <trigger> {string} </trigger>
+  <trigger> {property} </trigger>
 </pid>
 @endcode
 
+For the integration constant element, one can also supply the type attribute for
+what kind of integrator to be used, one of:
+
+- rect, for a rectangular integrator
+- trap, for a trapezoidal integrator
+- ab2, for a second order Adams Bashforth integrator
+- ab3, for a third order Adams Bashforth integrator
+
+For example,
+
+@code
+<pid name="fcs/heading-control">
+  <kp> 3 </kp>
+  <ki type="ab3"> 1 </ki>
+  <kd> 1 </kd>
+</pid>
+@endcode
+
+
+
 <h3>Configuration Parameters:</h3>
 <pre>
 
+  The values of kp, ki, and kd have slightly different interpretations depending on
+  whether the PID controller is a standard one, or an ideal/parallel one - with the latter
+  being the default.
+  
   kp      - Proportional constant, default value 0.
   ki      - Integrative constant, default value 0.
   kd      - Derivative constant, default value 0.
-  trigger - Property which is used to sense wind-up, optional.
+  trigger - Property which is used to sense wind-up, optional. Most often, the trigger
+            will be driven by the "saturated" property of a particular actuator. When
+            the relevant actuator has reached it's limits (if there are any, specified
+            by the <clipto> element) the automatically generated saturated property will
+            be greater than zero (true). If this property is used as the trigger for the
+            integrator, the integrator will not continue to integrate while the property
+            is still true (> 1), preventing wind-up.
+  pvdot   - The property to be used as the process variable time derivative. 
+
+
 
 </pre>
 
     @author Jon S. Berndt
-    @version $Revision: 1.12 $
+    @version $Revision: 1.13 $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -99,17 +132,26 @@ public:
   bool Run (void);
   void ResetPastStates(void) {Input_prev = Input_prev2 = Output = I_out_total = 0.0;}
 
+    /// These define the indices use to select the various integrators.
+  enum eIntegrateType {eNone = 0, eRectEuler, eTrapezoidal, eAdamsBashforth2, eAdamsBashforth3};
+
 private:
-  FGPropertyManager *Trigger;
   double Kp, Ki, Kd;
   double I_out_total;
   double Input_prev, Input_prev2;
   double KpPropertySign;
   double KiPropertySign;
   double KdPropertySign;
+
+  bool IsStandard;
+
+  eIntegrateType IntType;
+
+  FGPropertyManager *Trigger;
   FGPropertyManager* KpPropertyNode;
   FGPropertyManager* KiPropertyNode;
   FGPropertyManager* KdPropertyNode;
+  FGPropertyManager* ProcessVariableDot;
 
   void Debug(int from);
 };

src/FDM/JSBSim/models/flight_control/FGSensor.h

@@ -44,7 +44,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_SENSOR "$Id: FGSensor.h,v 1.20 2011/08/18 12:42:17 jberndt Exp $"
+#define ID_SENSOR "$Id: FGSensor.h,v 1.21 2012/01/08 12:39:14 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -74,6 +74,7 @@ Syntax:
     <max> number </max>
   </quantization>
   <drift_rate> number </drift_rate>
+  <gain> number </gain>
   <bias> number </bias>
   <delay> number < /delay>
 </sensor>
@@ -123,7 +124,7 @@ The delay element can specify a frame delay. The integer number provided is
 the number of frames to delay the output signal.
 
 @author Jon S. Berndt
-@version $Revision: 1.20 $
+@version $Revision: 1.21 $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

@@ -53,7 +53,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGEngine.cpp,v 1.48 2011/10/31 14:54:41 bcoconni Exp $";
+static const char *IdSrc = "$Id: FGEngine.cpp,v 1.50 2012/03/17 20:46:29 jentron Exp $";
 static const char *IdHdr = ID_ENGINE;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -248,10 +248,10 @@ bool FGEngine::LoadThruster(Element *thruster_element)
 
   thruster_filename = thruster_element->GetAttributeValue("file");
   if ( !thruster_filename.empty()) {
-    thruster_fullpathname = fullpath + thruster_filename + ".xml";
+    thruster_fullpathname = localpath + thruster_filename + ".xml";
     thruster_file.open(thruster_fullpathname.c_str());
     if ( !thruster_file.is_open()) {
-      thruster_fullpathname = localpath + thruster_filename + ".xml";
+      thruster_fullpathname = fullpath + thruster_filename + ".xml";
       thruster_file.open(thruster_fullpathname.c_str());
       if ( !thruster_file.is_open()) {
         cerr << "Could not open thruster file: " << thruster_filename << ".xml" << endl;

src/FDM/JSBSim/models/propulsion/FGEngine.h

@@ -55,7 +55,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_ENGINE "$Id: FGEngine.h,v 1.29 2011/11/10 12:06:14 jberndt Exp $"
+#define ID_ENGINE "$Id: FGEngine.h,v 1.35 2012/04/14 18:10:44 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -113,7 +113,7 @@ CLASS DOCUMENTATION
   documentation for engine and thruster classes.
 </pre>     
     @author Jon S. Berndt
-    @version $Id: FGEngine.h,v 1.29 2011/11/10 12:06:14 jberndt Exp $
+    @version $Id: FGEngine.h,v 1.35 2012/04/14 18:10:44 bcoconni Exp $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -149,6 +149,7 @@ public:
     vector <double> MixturePos;
     vector <double> PropAdvance;
     vector <bool> PropFeather;
+    FGColumnVector3 vXYZcg; // CG coordinates expressed in the structural frame
     double TotalDeltaT;
   };
 

src/FDM/JSBSim/models/propulsion/FGForce.h

@@ -66,7 +66,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_FORCE "$Id: FGForce.h,v 1.14 2011/10/31 14:54:41 bcoconni Exp $"
+#define ID_FORCE "$Id: FGForce.h,v 1.17 2012/04/01 17:05:51 bcoconni Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -215,7 +215,7 @@ and vMn, the moments, can be made directly. Otherwise, the usage is similar.<br>
 <br><br></p>
 
     @author Tony Peden
-    @version $Id: FGForce.h,v 1.14 2011/10/31 14:54:41 bcoconni Exp $
+    @version $Id: FGForce.h,v 1.17 2012/04/01 17:05:51 bcoconni Exp $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -321,14 +321,13 @@ protected:
   TransformType ttype;
   FGColumnVector3 vXYZn;
   FGColumnVector3 vActingXYZn;
+  FGMatrix33 mT;
 
 private:
   FGColumnVector3 vFb;
   FGColumnVector3 vM;
   FGColumnVector3 vDXYZ;
 
-  FGMatrix33 mT;
-
   void Debug(int from);
 };
 }

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

@@ -46,7 +46,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGNozzle.cpp,v 1.14 2011/08/03 03:21:06 jberndt Exp $";
+static const char *IdSrc = "$Id: FGNozzle.cpp,v 1.15 2012/03/18 15:48:35 jentron Exp $";
 static const char *IdHdr = ID_NOZZLE;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -90,7 +90,7 @@ double FGNozzle::Calculate(double vacThrust)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-string FGNozzle::GetThrusterLabels(int id, string delimeter)
+string FGNozzle::GetThrusterLabels(int id, const string& delimeter)
 {
   std::ostringstream buf;
 
@@ -101,7 +101,7 @@ string FGNozzle::GetThrusterLabels(int id, string delimeter)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-string FGNozzle::GetThrusterValues(int id, string delimeter)
+string FGNozzle::GetThrusterValues(int id, const string& delimeter)
 {
   std::ostringstream buf;
 

src/FDM/JSBSim/models/propulsion/FGNozzle.h

@@ -44,7 +44,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_NOZZLE "$Id: FGNozzle.h,v 1.9 2011/08/03 03:21:06 jberndt Exp $";
+#define ID_NOZZLE "$Id: FGNozzle.h,v 1.10 2012/03/18 15:48:36 jentron Exp $";
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -75,7 +75,7 @@ CLASS DOCUMENTATION
 
     All parameters MUST be specified.  
     @author Jon S. Berndt
-    @version $Id: FGNozzle.h,v 1.9 2011/08/03 03:21:06 jberndt Exp $
+    @version $Id: FGNozzle.h,v 1.10 2012/03/18 15:48:36 jentron Exp $
 */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -90,8 +90,8 @@ public:
   ~FGNozzle();
 
   double Calculate(double vacThrust);
-  string GetThrusterLabels(int id, string delimeter);
+  string GetThrusterLabels(int id, const string& delimeter);
-  string GetThrusterValues(int id, string delimeter);
+  string GetThrusterValues(int id, const string& delimeter);
 
 private:
 //  double PE;

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

@@ -50,7 +50,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGPiston.cpp,v 1.68 2011/10/11 15:13:34 jentron Exp $";
+static const char *IdSrc = "$Id: FGPiston.cpp,v 1.71 2012/04/07 01:50:54 jentron Exp $";
 static const char *IdHdr = ID_PISTON;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -85,6 +85,7 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
   MaxHP = 200;
   MinManifoldPressure_inHg = 6.5;
   MaxManifoldPressure_inHg = 28.5;
+  ManifoldPressureLag=1.0;
   ISFC = -1;
   volumetric_efficiency = 0.85;
   Bore = 5.125;
@@ -99,6 +100,9 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
   FMEPStatic = 46500;
   Cooling_Factor = 0.5144444;
   StaticFriction_HP = 1.5;
+  StarterGain = 1.;
+  StarterTorque = -1.;
+  StarterRPM = -1.;
 
   // These are internal program variables
 
@@ -120,6 +124,8 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
   bBoostOverride = false;
   bTakeoffBoost = false;
   TakeoffBoost = 0.0;   // Default to no extra takeoff-boost
+  BoostLossFactor = 0.0;   // Default to free boost
+  
   int i;
   for (i=0; i<FG_MAX_BOOST_SPEEDS; i++) {
     RatedBoost[i] = 0.0;
@@ -137,10 +143,12 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
 
   // Read inputs from engine data file where present.
 
-  if (el->FindElement("minmp")) // Should have ELSE statement telling default value used?
+  if (el->FindElement("minmp")) 
     MinManifoldPressure_inHg = el->FindElementValueAsNumberConvertTo("minmp","INHG");
   if (el->FindElement("maxmp"))
     MaxManifoldPressure_inHg = el->FindElementValueAsNumberConvertTo("maxmp","INHG");
+  if (el->FindElement("man-press-lag"))
+    ManifoldPressureLag = el->FindElementValueAsNumber("man-press-lag");
   if (el->FindElement("displacement"))
     Displacement = el->FindElementValueAsNumberConvertTo("displacement","IN3");
   if (el->FindElement("maxhp"))
@@ -179,6 +187,10 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
     Ram_Air_Factor  = el->FindElementValueAsNumber("ram-air-factor");
   if (el->FindElement("cooling-factor"))
     Cooling_Factor  = el->FindElementValueAsNumber("cooling-factor");
+  if (el->FindElement("starter-rpm"))
+    StarterRPM  = el->FindElementValueAsNumber("starter-rpm");
+  if (el->FindElement("starter-torque"))
+    StarterTorque  = el->FindElementValueAsNumber("starter-torque");
   if (el->FindElement("dynamic-fmep"))
     FMEPDynamic= el->FindElementValueAsNumberConvertTo("dynamic-fmep","PA");
   if (el->FindElement("static-fmep"))
@@ -193,6 +205,8 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
       BoostManual = (int)el->FindElementValueAsNumber("boostmanual");
     if (el->FindElement("takeoffboost"))
       TakeoffBoost = el->FindElementValueAsNumberConvertTo("takeoffboost", "PSI");
+    if (el->FindElement("boost-loss-factor"))
+      BoostLossFactor = el->FindElementValueAsNumber("boost-loss-factor");
     if (el->FindElement("ratedboost1"))
       RatedBoost[0] = el->FindElementValueAsNumberConvertTo("ratedboost1", "PSI");
     if (el->FindElement("ratedboost2"))
@@ -234,7 +248,13 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
     }
   }
 
-  StarterHP = sqrt(MaxHP) * 0.4;
+
+  volumetric_efficiency_reduced = volumetric_efficiency;
+
+  if(StarterRPM < 0.) StarterRPM = 2*IdleRPM;
+  if(StarterTorque < 0)
+    StarterTorque = (MaxHP)*0.4; //just a wag.
+
   displacement_SI = Displacement * in3tom3;
   RatedMeanPistonSpeed_fps =  ( MaxRPM * Stroke) / (360); // AKA 2 * (RPM/60) * ( Stroke / 12) or 2NS
 
@@ -318,8 +338,12 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
   base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNumber);
   property_name = base_property_name + "/power-hp";
   PropertyManager->Tie(property_name, &HP);
+  property_name = base_property_name + "/friction-hp";
+  PropertyManager->Tie(property_name, &StaticFriction_HP);
   property_name = base_property_name + "/bsfc-lbs_hphr";
   PropertyManager->Tie(property_name, &ISFC);
+  property_name = base_property_name + "/starter-norm";
+  PropertyManager->Tie(property_name, &StarterGain);
   property_name = base_property_name + "/volumetric-efficiency";
   PropertyManager->Tie(property_name, &volumetric_efficiency);
   property_name = base_property_name + "/map-pa";
@@ -342,6 +366,12 @@ FGPiston::FGPiston(FGFDMExec* exec, Element* el, int engine_number, struct Input
   PropertyManager->Tie(property_name, this, &FGPiston::getOilPressure_psi);
   property_name = base_property_name + "/egt-degF";
   PropertyManager->Tie(property_name, this, &FGPiston::getExhaustGasTemp_degF);
+  if(BoostLossFactor > 0.0) {
+    property_name = base_property_name + "/boostloss-factor";
+    PropertyManager->Tie(property_name, &BoostLossFactor);
+    property_name = base_property_name + "/boostloss-hp";
+    PropertyManager->Tie(property_name, &BoostLossHP);
+  }
 
   // Set up and sanity-check the turbo/supercharging configuration based on the input values.
   if (TakeoffBoost > RatedBoost[0]) bTakeoffBoost = true;
@@ -419,6 +449,7 @@ void FGPiston::ResetToIC(void)
   Thruster->SetRPM(0.0);
   RPM = 0.0;
   OilPressure_psi = 0.0;
+  BoostLossHP = 0.;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -490,8 +521,8 @@ double FGPiston::CalcFuelNeed(void)
 int FGPiston::InitRunning(void)
 {
   Magnetos=3;
-  in.MixtureCmd[EngineNumber] = in.PressureRatio/1.3;
+  in.MixtureCmd[EngineNumber] = in.PressureRatio*1.3;
-  in.MixturePos[EngineNumber] = in.PressureRatio/1.3;
+  in.MixturePos[EngineNumber] = in.PressureRatio*1.3;
   Thruster->SetRPM( 2.0*IdleRPM/Thruster->GetGearRatio() );
   Running = true;
   return 1;
@@ -525,42 +556,30 @@ void FGPiston::doEngineStartup(void)
   if ((Magnetos == 1) || (Magnetos > 2)) Magneto_Left = true;
   if (Magnetos > 1)  Magneto_Right = true;
 
-  // Assume we have fuel for now
+// We will 'run' with any fuel flow. If there is not enough fuel to make power it will show in doEnginePower
-  fuel = !Starved;
+  fuel = FuelFlowRate > 0.0 ? 1 : 0;
 
   // Check if we are turning the starter motor
   if (Cranking != Starter) {
     // This check saves .../cranking from getting updated every loop - they
     // only update when changed.
     Cranking = Starter;
-    crank_counter = 0;
   }
 
-  if (Cranking) crank_counter++;  //Check mode of engine operation
 
-  if (!Running && spark && fuel) {  // start the engine if revs high enough
+  // Cut the engine *power* - Note: the engine will continue to
-    if (Cranking) {
+  // spin depending on prop Ixx and freestream velocity
-      if ((RPM > IdleRPM*0.8) && (crank_counter > 175)) // Add a little delay to startup
+
-        Running = true;                         // on the starter
+  if ( Running ) {
-    } else {
+    if (!spark || !fuel)    Running = false;
-      if (RPM > IdleRPM*0.8)                            // This allows us to in-air start
+    if (RPM < IdleRPM*0.8 ) Running = false;
-        Running = true;                         // when windmilling
+  } else { // !Running  
+    if ( spark && fuel) {     // start the engine if revs high enough
+      if (RPM > IdleRPM*0.8)  // This allows us to in-air start
+        Running = true;       // when windmilling
     }
   }
 
-  // Cut the engine *power* - Note: the engine may continue to
-  // spin if the prop is in a moving airstream
-
-  if ( Running && (!spark || !fuel) ) Running = false;
-
-  // Check for stalling (RPM = 0).
-  if (Running) {
-    if (RPM == 0) {
-      Running = false;
-    } else if ((RPM <= IdleRPM *0.8 ) && (Cranking)) {
-      Running = false;
-    }
-  }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -578,7 +597,7 @@ void FGPiston::doEngineStartup(void)
 
 void FGPiston::doBoostControl(void)
 {
-  if(BoostManual) {
+  if(bBoostManual) {
     if(BoostSpeed > BoostSpeeds-1) BoostSpeed = BoostSpeeds-1;
     if(BoostSpeed < 0) BoostSpeed = 0;
   } else {
@@ -608,7 +627,7 @@ void FGPiston::doBoostControl(void)
  * Inputs: p_amb, Throttle,
  *         MeanPistonSpeed_fps, dt
  *
- * Outputs: MAP, ManifoldPressure_inHg, TMAP
+ * Outputs: MAP, ManifoldPressure_inHg, TMAP, BoostLossHP
  */
 
 void FGPiston::doMAP(void)
@@ -618,9 +637,9 @@ void FGPiston::doMAP(void)
 
   double map_coefficient = Ze/(Ze+Z_airbox+Zt);
 
-  // Add a one second lag to manifold pressure changes
+  // Add a variable lag to manifold pressure changes
-  double dMAP=0;
+  double dMAP=(TMAP - p_ram * map_coefficient);
-  dMAP = (TMAP - p_ram * map_coefficient) * TotalDeltaT;
+  if (ManifoldPressureLag > TotalDeltaT) dMAP *= TotalDeltaT/ManifoldPressureLag;
 
   TMAP -=dMAP;
 
@@ -646,15 +665,26 @@ void FGPiston::doMAP(void)
     double boost_factor = (( BoostMul[BoostSpeed] - 1 ) / RatedRPM[BoostSpeed] ) * RPM + 1;
     MAP = TMAP * boost_factor;
     // Now clip the manifold pressure to BCV or Wastegate setting.
-    if (bTakeoffPos) {
+    if(!bBoostOverride) {
-      if (MAP > TakeoffMAP[BoostSpeed]) MAP = TakeoffMAP[BoostSpeed];
+      if (bTakeoffPos) {
-    } else {
+        if (MAP > TakeoffMAP[BoostSpeed]) MAP = TakeoffMAP[BoostSpeed];
-      if (MAP > RatedMAP[BoostSpeed]) MAP = RatedMAP[BoostSpeed];
+      } else {
+        if (MAP > RatedMAP[BoostSpeed]) MAP = RatedMAP[BoostSpeed];
+      }
     }
   } else {
       MAP = TMAP;
   }
 
+  if( BoostLossFactor > 0.0 )
+  {
+      double gamma = 1.414; // specific heat constants
+      double Nstage = 1; // Nstage is the number of boost stages.
+      BoostLossHP = ((Nstage * TMAP * v_dot_air * gamma) / (gamma - 1)) * (pow((MAP/TMAP),((gamma-1)/(Nstage * gamma))) - 1) * BoostLossFactor / 745.7 ; // 745.7 convert watt to hp
+  } else {
+      BoostLossHP = 0;
+  }
+  
   // And set the value in American units as well
   ManifoldPressure_inHg = MAP / inhgtopa;
 }
@@ -670,7 +700,7 @@ void FGPiston::doMAP(void)
  *
  * TODO: Model inlet manifold air temperature.
  *
- * Outputs: rho_air, m_dot_air
+ * Outputs: rho_air, m_dot_air, volumetric_efficiency_reduced
  */
 
 void FGPiston::doAirFlow(void)
@@ -678,11 +708,14 @@ void FGPiston::doAirFlow(void)
   double gamma = 1.3; // specific heat constants
 // loss of volumentric efficiency due to difference between MAP and exhaust pressure
 // Eq 6-10 from The Internal Combustion Engine - Charles Taylor Vol 1
-  double ve =((gamma-1)/gamma) +( CompressionRatio -(p_amb/MAP))/(gamma*( CompressionRatio - 1));
+  double mratio = MAP < 1. ? CompressionRatio : p_amb/MAP;
+  if (mratio > CompressionRatio) mratio = CompressionRatio;
+  double ve =((gamma-1)/gamma) +( CompressionRatio -(mratio))/(gamma*( CompressionRatio - 1));
 
   rho_air = p_amb / (R_air * T_amb);
   double swept_volume = (displacement_SI * (RPM/60)) / 2;
-  double v_dot_air = swept_volume * volumetric_efficiency *ve;
+  volumetric_efficiency_reduced = volumetric_efficiency *ve;
+  v_dot_air = swept_volume * volumetric_efficiency_reduced;
 
   double rho_air_manifold = MAP / (R_air * T_amb);
   m_dot_air = v_dot_air * rho_air_manifold;
@@ -729,14 +762,12 @@ void FGPiston::doFuelFlow(void)
 
 void FGPiston::doEnginePower(void)
 {
-  IndicatedHorsePower = 0;
+  IndicatedHorsePower = -StaticFriction_HP;
   FMEP = 0;
   if (Running) {
-    // FIXME: this needs to be generalized
+    double ME, power;  // Convienience term for use in the calculations
-    double ME, percent_RPM, power;  // Convienience term for use in the calculations
     ME = Mixture_Efficiency_Correlation->GetValue(m_dot_fuel/m_dot_air);
 
-    percent_RPM = RPM/MaxRPM;
 // Guestimate engine friction losses from Figure 4.4 of "Engines: An Introduction", John Lumley
     FMEP = (-FMEPDynamic * MeanPistonSpeed_fps * fttom - FMEPStatic);
 
@@ -745,27 +776,26 @@ void FGPiston::doEnginePower(void)
     if ( Magnetos != 3 ) power *= SparkFailDrop;
 
 
-    IndicatedHorsePower = (FuelFlow_pph / ISFC )* ME * power;
+    IndicatedHorsePower = (FuelFlow_pph / ISFC )* ME * power - StaticFriction_HP; //FIXME static friction should depend on oil temp and configuration;
 
   } else {
     // Power output when the engine is not running
+    double torque, k_torque, rpm;  // Convienience term for use in the calculations
+    
+    rpm = RPM < 1.0 ? 1.0 : RPM;
     if (Cranking) {
-      if (RPM < 10) {
+      if(RPM<StarterRPM) k_torque = 1.0-RPM/(StarterRPM);
-        IndicatedHorsePower = StarterHP;
+      else k_torque = 0;
-      } else if (RPM < IdleRPM*0.8) {
+      torque = StarterTorque*k_torque*StarterGain;
-        IndicatedHorsePower = StarterHP + ((IdleRPM*0.8 - RPM) / 8.0);
+      IndicatedHorsePower = torque * rpm / 5252;
-        // This is a guess - would be nice to find a proper starter moter torque curve
+     } 
-      } else {
-        IndicatedHorsePower = StarterHP;
-      }
-    }
   }
 
   // Constant is (1/2) * 60 * 745.7
   // (1/2) convert cycles, 60 minutes to seconds, 745.7 watts to hp.
   double pumping_hp = ((PMEP + FMEP) * displacement_SI * RPM)/(Cycles*22371);
 
-  HP = IndicatedHorsePower + pumping_hp - StaticFriction_HP; //FIXME static friction should depend on oil temp and configuration
+HP = IndicatedHorsePower + pumping_hp - BoostLossHP;
 //  cout << "pumping_hp " <<pumping_hp << FMEP << PMEP <<endl;
   PctPower = HP / MaxHP ;
 //  cout << "Power = " << HP << "  RPM = " << RPM << "  Running = " << Running << "  Cranking = " << Cranking << endl;
@@ -989,7 +1019,7 @@ void FGPiston::Debug(int from)
       cout << "      Bore: "                << Bore                     << endl;
       cout << "      Stroke: "              << Stroke                   << endl;
       cout << "      Cylinders: "           << Cylinders                << endl;
-      cout << "      Cylinders Head Mass: " <<CylinderHeadMass          << endl;
+      cout << "      Cylinders Head Mass: " << CylinderHeadMass         << endl;
       cout << "      Compression Ratio: "   << CompressionRatio         << endl;
       cout << "      MaxHP: "               << MaxHP                    << endl;
       cout << "      Cycles: "              << Cycles                   << endl;
@@ -1003,6 +1033,9 @@ void FGPiston::Debug(int from)
       cout << "      Dynamic FMEP Factor: " << FMEPDynamic << endl;
       cout << "      Static FMEP Factor: " << FMEPStatic << endl;
 
+      cout << "      Starter Motor Torque: " << StarterTorque << endl;
+      cout << "      Starter Motor RPM:    " << StarterRPM << endl;
+
       cout << endl;
       cout << "      Combustion Efficiency table:" << endl;
       Lookup_Combustion_Efficiency->Print();

src/FDM/JSBSim/models/propulsion/FGPiston.h

@@ -46,7 +46,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_PISTON "$Id: FGPiston.h,v 1.32 2011/10/11 16:16:16 jentron Exp $";
+#define ID_PISTON "$Id: FGPiston.h,v 1.35 2012/04/07 01:50:54 jentron Exp $";
 #define FG_MAX_BOOST_SPEEDS 3
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -81,6 +81,9 @@ CLASS DOCUMENTATION
   <air-intake-impedance-factor> {number} </air-intake-impedance-factor>
   <ram-air-factor> {number} </ram-air-factor>
   <cooling-factor> {number} </cooling-factor>
+  <man-press-lag> {number} </man-press-lag>
+  <starter-torque> {number} </starter-torque> 
+  <starter-rpm> {number} </starter-rpm> 
   <cylinder-head-mass unit="{KG | LBS}"> {number} </cylinder-head-mass>
   <bsfc unit="{LBS/HP*HR | "KG/KW*HR"}"> {number} </bsfc>
   <volumetric-efficiency> {number} </volumetric-efficiency>
@@ -89,6 +92,7 @@ CLASS DOCUMENTATION
   <numboostspeeds> {number} </numboostspeeds>
   <boostoverride> {0 | 1} </boostoverride>
   <boostmanual> {0 | 1} </boostmanual>
+  <boost-loss-factor> {number} </boost-loss-factor>
   <ratedboost1 unit="{INHG | PA | ATM}"> {number} </ratedboost1>
   <ratedpower1 unit="{HP | WATTS}"> {number} </ratedpower1>
   <ratedrpm1> {number} </ratedrpm1>
@@ -112,6 +116,12 @@ Basic parameters:
 - \b maxmp - this value is the nominal maximum manifold pressure at sea-level
       without boost. Along with maxrpm it determines the resistance of the
       aircraft intake system. Overridden by air-intake-impedance-factor
+- \b man-press-lag - Delay in seconds for manifold pressure changes to take effect
+- \b starter-torque - A value specifing the zero RPM torque in lb*ft the starter motor
+      provides. Current default value is 40% of the horse power value.
+- \b starter-rpm - A value specifing the maximum RPM the unloaded starter motor
+      can achieve. Loads placed on the engine by the propeller and throttle will
+      further limit RPM achieved in practice.
 - \b idlerpm - this value affects the throttle fall off and the engine stops
       running if it is slowed below 80% of this value. The engine starts
       running when it reaches 80% of this value.
@@ -163,7 +173,10 @@ Supercharge parameters:
       supercharger speeds.  Merlin XII had 1 speed, Merlin 61 had 2, a late
       Griffon engine apparently had 3.  No known engine more than 3, although
       some German engines had continuously variable-speed superchargers.
-- \b boostoverride - unused
+- \b boostoverride - whether or not to clip output to the wastegate value
+- \b boost-loss-factor - zero (or not present) for 'free' supercharging. A value entered
+      will be used as a multiplier to the power required to compress the input air. Typical 
+      value should be 1.15 to 1.20.
 - \b boostmanual - whether a multispeed supercharger will manually or
       automatically shift boost speeds. On manual shifting the boost speeds is
       accomplished by controlling the property propulsion/engine/boostspeed.
@@ -198,7 +211,7 @@ boostspeed they refer to:
     @author David Megginson (initial porting and additional code)
     @author Ron Jensen (additional engine code)
     @see Taylor, Charles Fayette, "The Internal Combustion Engine in Theory and Practice"
-    @version $Id: FGPiston.h,v 1.32 2011/10/11 16:16:16 jentron Exp $
+    @version $Id: FGPiston.h,v 1.35 2012/04/07 01:50:54 jentron Exp $
   */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -239,10 +252,12 @@ private:
   int crank_counter;
 
   double IndicatedHorsePower;
+  double IndicatedPower;
   double PMEP;
   double FMEP;
   double FMEPDynamic;
   double FMEPStatic;
+  double T_Intake;
 
   void doEngineStartup(void);
   void doBoostControl(void);
@@ -279,6 +294,7 @@ private:
   double MinManifoldPressure_inHg; // Inches Hg
   double MaxManifoldPressure_inHg; // Inches Hg
   double MaxManifoldPressure_Percent; // MaxManifoldPressure / 29.92
+  double ManifoldPressureLag;      // Manifold Pressure delay in seconds.
   double Displacement;             // cubic inches
   double displacement_SI;          // cubic meters
   double MaxHP;                    // horsepower
@@ -298,7 +314,9 @@ private:
   double RatedMeanPistonSpeed_fps; // ft/sec derived from MaxRPM and stroke.
   double Ram_Air_Factor;           // number
 
-  double StarterHP;                // initial horsepower of starter motor
+  double StarterTorque;// Peak Torque of the starter motor
+  double StarterRPM;   // Peak RPM of the starter motor
+  double StarterGain;  // control the torque of the starter motor.
   int BoostSpeeds;  // Number of super/turbocharger boost speeds - zero implies no turbo/supercharging.
   int BoostSpeed;   // The current boost-speed (zero-based).
   bool Boosted;     // Set true for boosted engine.
@@ -322,6 +340,7 @@ private:
   double RatedMAP[FG_MAX_BOOST_SPEEDS]; // Rated manifold absolute pressure [Pa] (BCV clamp)
   double TakeoffMAP[FG_MAX_BOOST_SPEEDS];   // Takeoff setting manifold absolute pressure [Pa] (BCV clamp)
   double BoostSwitchHysteresis; // Pa.
+  double BoostLossFactor; // multiplier for HP consumed by the supercharger
 
   double minMAP;  // Pa
   double maxMAP;  // Pa
@@ -348,11 +367,14 @@ private:
   //
   double rho_air;
   double volumetric_efficiency;
+  double volumetric_efficiency_reduced;
   double map_coefficient;
   double m_dot_air;
+  double v_dot_air;
   double equivalence_ratio;
   double m_dot_fuel;
   double HP;
+  double BoostLossHP;
   double combustion_efficiency;
   double ExhaustGasTemp_degK;
   double EGT_degC;

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

@@ -45,7 +45,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGPropeller.cpp,v 1.41 2011/11/17 21:07:30 jentron Exp $";
+static const char *IdSrc = "$Id: FGPropeller.cpp,v 1.44 2012/04/29 13:10:46 bcoconni Exp $";
 static const char *IdHdr = ID_PROPELLER;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -228,16 +228,26 @@ double FGPropeller::Calculate(double EnginePower)
   // Induced velocity in the propeller disk area. This formula is obtained
   // from momentum theory - see B. W. McCormick, "Aerodynamics, Aeronautics,
   // and Flight Mechanics" 1st edition, eqn. 6.15 (propeller analysis chapter).
-  // Vinduced = 0.5 * (-Vel + sqrt(Vel*Vel + 2.0*Thrust/(rho*Area)))
   // Since Thrust and Vel can both be negative we need to adjust this formula
   // To handle sign (direction) separately from magnitude.
   double Vel2sum = Vel*abs(Vel) + 2.0*Thrust/(rho*Area);
-
+  
   if( Vel2sum > 0.0)
     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)
+    Vinduced = 0.0; // We cannot calculate the induced velocity so let's assume it is zero.
+    
   // 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
   // and the relative wind that goes through the propeller disk.

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

@@ -47,7 +47,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGRocket.cpp,v 1.26 2011/08/04 13:45:42 jberndt Exp $";
+static const char *IdSrc = "$Id: FGRocket.cpp,v 1.27 2012/04/08 15:19:08 jberndt Exp $";
 static const char *IdHdr = ID_ROCKET;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -55,7 +55,7 @@ CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 FGRocket::FGRocket(FGFDMExec* exec, Element *el, int engine_number, struct Inputs& input)
-  : FGEngine(exec, el, engine_number, input)
+  : FGEngine(exec, el, engine_number, input), isp_function(0L)
 {
   Type = etRocket;
   Element* thrust_table_element = 0;
@@ -67,7 +67,8 @@ FGRocket::FGRocket(FGFDMExec* exec, Element *el, int engine_number, struct Input
   TotalPropellantExpended = 0.0;
   FuelFlowRate = FuelExpended = 0.0;
   OxidizerFlowRate = OxidizerExpended = 0.0;
-  SLOxiFlowMax = SLFuelFlowMax = 0.0;
+  SLOxiFlowMax = SLFuelFlowMax = PropFlowMax = 0.0;
+  MxR = 0.0;
   BuildupTime = 0.0;
   It = 0.0;
   ThrustVariation = 0.0;
@@ -78,19 +79,51 @@ FGRocket::FGRocket(FGFDMExec* exec, Element *el, int engine_number, struct Input
    MinThrottle = 0.0;
    MaxThrottle = 1.0;
 
-  if (el->FindElement("isp"))
+  string base_property_name = CreateIndexedPropertyName("propulsion/engine", EngineNumber);
+
+  std::stringstream strEngineNumber;
+  strEngineNumber << EngineNumber;
+
+  Element* isp_el = el->FindElement("isp");
+  Element* isp_func_el=0;
+
+  bindmodel(); // Bind model properties first, since they might be needed in functions.
+
+  // Specific impulse may be specified as a constant value or as a function - perhaps as a function of mixture ratio.
+  if (isp_el) {
+    isp_func_el = isp_el->FindElement("function");
+    if (isp_func_el) {
+      isp_function = new FGFunction(exec->GetPropertyManager(),isp_func_el, strEngineNumber.str());
+    } else {
     Isp = el->FindElementValueAsNumber("isp");
+    }
+  } else {
+    throw("Specific Impulse <isp> must be specified for a rocket engine");
+  }
+  
   if (el->FindElement("builduptime"))
     BuildupTime = el->FindElementValueAsNumber("builduptime");
   if (el->FindElement("maxthrottle"))
     MaxThrottle = el->FindElementValueAsNumber("maxthrottle");
   if (el->FindElement("minthrottle"))
     MinThrottle = el->FindElementValueAsNumber("minthrottle");
-  if (el->FindElement("slfuelflowmax"))
-    SLFuelFlowMax = el->FindElementValueAsNumberConvertTo("slfuelflowmax", "LBS/SEC");
-  if (el->FindElement("sloxiflowmax"))
-    SLOxiFlowMax = el->FindElementValueAsNumberConvertTo("sloxiflowmax", "LBS/SEC");
 
+  if (el->FindElement("slfuelflowmax")) {
+    SLFuelFlowMax = el->FindElementValueAsNumberConvertTo("slfuelflowmax", "LBS/SEC");
+    if (el->FindElement("sloxiflowmax")) {
+    SLOxiFlowMax = el->FindElementValueAsNumberConvertTo("sloxiflowmax", "LBS/SEC");
+    }
+    PropFlowMax = SLOxiFlowMax + SLFuelFlowMax;
+    MxR = SLOxiFlowMax/SLFuelFlowMax;
+  } else if (el->FindElement("propflowmax")) {
+    PropFlowMax = el->FindElementValueAsNumberConvertTo("propflowmax", "LBS/SEC");
+    // Mixture ratio may be specified here, but it can also be specified as a function or via property
+    if (el->FindElement("mixtureratio")) {
+      MxR = el->FindElementValueAsNumber("mixtureratio");
+    }
+  }
+
+  if (isp_function) Isp = isp_function->GetValue(); // cause Isp function to be executed if present.
   // If there is a thrust table element, this is a solid propellant engine.
   thrust_table_element = el->FindElement("thrust_table");
   if (thrust_table_element) {
@@ -106,7 +139,6 @@ FGRocket::FGRocket(FGFDMExec* exec, Element *el, int engine_number, struct Input
     }
   }
 
-  bindmodel();
 
   Debug(0);
 }
@@ -129,6 +161,9 @@ void FGRocket::Calculate(void)
 
   PropellantFlowRate = (FuelExpended + OxidizerExpended)/in.TotalDeltaT;
   TotalPropellantExpended += FuelExpended + OxidizerExpended;
+  // If Isp has been specified as a function, override the value of Isp to that, otherwise
+  // assume a constant value is given.
+  if (isp_function) Isp = isp_function->GetValue();
 
   // If there is a thrust table, it is a function of propellant burned. The
   // engine is started when the throttle is advanced to 1.0. After that, it
@@ -189,7 +224,8 @@ double FGRocket::CalcFuelNeed(void)
     FuelFlowRate = VacThrust/Isp;   // This calculates wdot (weight flow rate in lbs/sec)
     FuelFlowRate /= (1 + TotalIspVariation);
   } else {
-    FuelFlowRate = SLFuelFlowMax*PctPower;
+    SLFuelFlowMax = PropFlowMax / (1 + MxR);
+    FuelFlowRate = SLFuelFlowMax * PctPower;
   }
 
   FuelExpended = FuelFlowRate * in.TotalDeltaT; // For this time step ...
@@ -200,6 +236,7 @@ double FGRocket::CalcFuelNeed(void)
 
 double FGRocket::CalcOxidizerNeed(void)
 {
+  SLOxiFlowMax = PropFlowMax * MxR / (1 + MxR);
   OxidizerFlowRate = SLOxiFlowMax * PctPower;
   OxidizerExpended = OxidizerFlowRate * in.TotalDeltaT;
   return OxidizerExpended;
@@ -252,6 +289,12 @@ void FGRocket::bindmodel()
   } else { // Liquid rocket motor
     property_name = base_property_name + "/oxi-flow-rate-pps";
     PropertyManager->Tie( property_name.c_str(), this, &FGRocket::GetOxiFlowRate);
+    property_name = base_property_name + "/mixture-ratio";
+    PropertyManager->Tie( property_name.c_str(), this, &FGRocket::GetMixtureRatio,
+                                                       &FGRocket::SetMixtureRatio);
+    property_name = base_property_name + "/isp";
+    PropertyManager->Tie( property_name.c_str(), this, &FGRocket::GetIsp,
+                                                       &FGRocket::SetIsp);
   }
 }
 

src/FDM/JSBSim/models/propulsion/FGRocket.h

@@ -40,13 +40,14 @@ INCLUDES
 
 #include "FGEngine.h"
 #include "math/FGTable.h"
+#include "math/FGFunction.h"
 #include "input_output/FGXMLElement.h"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_ROCKET "$Id: FGRocket.h,v 1.17 2011/08/04 13:45:42 jberndt Exp $"
+#define ID_ROCKET "$Id: FGRocket.h,v 1.18 2012/04/08 15:19:08 jberndt Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -118,7 +119,7 @@ for the rocket engine to be throttle up to 1. At that time, the solid rocket
 fuel begins burning and thrust is provided.
 
     @author Jon S. Berndt
-    $Id: FGRocket.h,v 1.17 2011/08/04 13:45:42 jberndt Exp $
+    $Id: FGRocket.h,v 1.18 2012/04/08 15:19:08 jberndt Exp $
     @see FGNozzle,
     FGThruster,
     FGForce,
@@ -168,6 +169,14 @@ public:
 
   double GetOxiFlowRate(void) const {return OxidizerFlowRate;}
 
+  double GetMixtureRatio(void) const {return MxR;}
+
+  double GetIsp(void) const {return Isp;}
+
+  void SetMixtureRatio(double mix) {MxR = mix;}
+
+  void SetIsp(double isp) {Isp = isp;}
+
   std::string GetEngineLabels(const std::string& delimiter);
   std::string GetEngineValues(const std::string& delimiter);
 
@@ -216,11 +225,13 @@ private:
   double OxidizerExpended;
   double TotalPropellantExpended;
   double SLOxiFlowMax;
+  double PropFlowMax;
   double OxidizerFlowRate;
   double PropellantFlowRate;
   bool Flameout;
   double BuildupTime;
   FGTable* ThrustTable;
+  FGFunction* isp_function;
 
   void Debug(int from);
 };

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

@@ -36,27 +36,27 @@ HISTORY
 01/10/11  T.Kreitler changed to single rotor model
 03/06/11  T.Kreitler added brake, clutch, and experimental free-wheeling-unit,
                      reasonable estimate for inflowlag
+02/05/12  T.Kreitler brake, clutch, and FWU now in FGTransmission, 
+                     downwash angles relate to shaft orientation
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#include <iostream>
 #include <sstream>
 
 #include "FGRotor.h"
-#include "input_output/FGXMLElement.h"
 #include "models/FGMassBalance.h"
 #include "models/FGPropulsion.h" // to get the GearRatio from a linked rotor
 
 using std::cerr;
+using std::cout;
 using std::endl;
 using std::ostringstream;
-using std::cout;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGRotor.cpp,v 1.18 2011/10/15 21:30:28 jentron Exp $";
+static const char *IdSrc = "$Id: FGRotor.cpp,v 1.20 2012/03/18 15:48:36 jentron Exp $";
 static const char *IdHdr = ID_ROTOR;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -69,116 +69,6 @@ static inline double sqr(double x) { return x*x; }
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-// Note: The FGTransmission class is currently carried 'pick-a-pack' by
-// FGRotor.
-
-FGTransmission::FGTransmission(FGFDMExec *exec, int num) :
-  FreeWheelTransmission(1.0),
-  ThrusterMoment(1.0), EngineMoment(1.0), EngineFriction(0.0),
-  ClutchCtrlNorm(1.0), BrakeCtrlNorm(0.0), MaxBrakePower(0.0),
-  EngineRPM(0.0), ThrusterRPM(0.0)
-{
-  double dt;
-  PropertyManager = exec->GetPropertyManager();
-  dt = exec->GetDeltaT();
-
-  // avoid too abrupt changes in transmission
-  FreeWheelLag = Filter(200.0,dt);
-  BindModel(num);
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-FGTransmission::~FGTransmission(){
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-// basically P = Q*w and Q_Engine + (-Q_Rotor) = J * dw/dt, J = Moment
-//
-void FGTransmission::Calculate(double EnginePower, double ThrusterTorque, double dt) {
-
-  double coupling = 1.0, coupling_sq = 1.0;
-  double fw_mult = 1.0;
-
-  double d_omega = 0.0, engine_d_omega = 0.0, thruster_d_omega = 0.0; // relative changes
-
-  double engine_omega = rpm_to_omega(EngineRPM);
-  double safe_engine_omega = engine_omega < 1e-1 ? 1e-1 : engine_omega;
-  double engine_torque = EnginePower / safe_engine_omega;
-
-  double thruster_omega = rpm_to_omega(ThrusterRPM);
-  double safe_thruster_omega = thruster_omega < 1e-1 ? 1e-1 : thruster_omega;
-
-  engine_torque  -= EngineFriction / safe_engine_omega;
-  ThrusterTorque += Constrain(0.0, BrakeCtrlNorm, 1.0) * MaxBrakePower / safe_thruster_omega;
-
-  // would the FWU release ?
-  engine_d_omega = engine_torque/EngineMoment * dt;
-  thruster_d_omega =  - ThrusterTorque/ThrusterMoment * dt;
-
-  if ( thruster_omega+thruster_d_omega > engine_omega+engine_d_omega ) {
-    // don't drive the engine
-    FreeWheelTransmission = 0.0;
-  } else {
-    FreeWheelTransmission = 1.0;
-  }
-
-  fw_mult = FreeWheelLag.execute(FreeWheelTransmission);
-  coupling = fw_mult * Constrain(0.0, ClutchCtrlNorm, 1.0);
-
-  if (coupling < 0.999999) { // are the separate calculations needed ?
-    // assume linear transfer 
-    engine_d_omega   =
-       (engine_torque - ThrusterTorque*coupling)/(ThrusterMoment*coupling + EngineMoment) * dt;
-    thruster_d_omega = 
-       (engine_torque*coupling - ThrusterTorque)/(ThrusterMoment + EngineMoment*coupling) * dt;
-
-    EngineRPM += omega_to_rpm(engine_d_omega);
-    ThrusterRPM += omega_to_rpm(thruster_d_omega);
-
-    // simulate friction in clutch
-    coupling_sq = coupling*coupling;
-    EngineRPM   = (1.0-coupling_sq) * EngineRPM    + coupling_sq * 0.02 * (49.0*EngineRPM + ThrusterRPM);
-    ThrusterRPM = (1.0-coupling_sq) * ThrusterRPM  + coupling_sq * 0.02 * (EngineRPM + 49.0*ThrusterRPM);
-
-    // enforce equal rpm
-    if ( fabs(EngineRPM-ThrusterRPM) < 1e-3 ) {
-      EngineRPM = ThrusterRPM = 0.5 * (EngineRPM+ThrusterRPM);
-    }
-  } else {
-    d_omega = (engine_torque - ThrusterTorque)/(ThrusterMoment + EngineMoment) * dt;
-    EngineRPM = ThrusterRPM += omega_to_rpm(d_omega);
-  }
-
-  // nothing will turn backward
-  if (EngineRPM < 0.0 ) EngineRPM = 0.0;
-  if (ThrusterRPM < 0.0 ) ThrusterRPM = 0.0;
-
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-bool FGTransmission::BindModel(int num)
-{
-  string property_name, base_property_name;
-  base_property_name = CreateIndexedPropertyName("propulsion/engine", num);
-
-  property_name = base_property_name + "/brake-ctrl-norm";
-  PropertyManager->Tie( property_name.c_str(), this, &FGTransmission::GetBrakeCtrl, &FGTransmission::SetBrakeCtrl);
-  property_name = base_property_name + "/free-wheel-transmission";
-  PropertyManager->Tie( property_name.c_str(), this, &FGTransmission::GetFreeWheelTransmission);
-
-  return true;
-}
-
-
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 // Constructor
 
 FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
@@ -190,7 +80,7 @@ FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
     BladeChord(0.0), LiftCurveSlope(0.0), BladeTwist(0.0), HingeOffset(0.0),
     BladeFlappingMoment(0.0), BladeMassMoment(0.0), PolarMoment(0.0),
     InflowLag(0.0), TipLossB(0.0),
-    GroundEffectExp(0.0), GroundEffectShift(0.0),
+    GroundEffectExp(0.0), GroundEffectShift(0.0), GroundEffectScaleNorm(1.0),
     LockNumberByRho(0.0), Solidity(0.0),            // derived parameters
     RPM(0.0), Omega(0.0),                           // dynamic values
     beta_orient(0.0),
@@ -206,6 +96,7 @@ FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
 {
   FGColumnVector3 location(0.0, 0.0, 0.0), orientation(0.0, 0.0, 0.0);
   Element *thruster_element;
+  double engine_power_est = 0.0;
 
   // initialise/set remaining variables
   SetTransformType(FGForce::tCustom);
@@ -246,7 +137,7 @@ FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
 
   SetLocation(location);
   SetAnglesToBody(orientation);
-  InvTransform = Transform().Transposed();
+  InvTransform = Transform().Transposed(); // body to custom/native
 
   // wire controls
   ControlMap = eMainCtrl;
@@ -283,41 +174,31 @@ FGRotor::FGRotor(FGFDMExec *exec, Element* rotor_element, int num)
     }
   }
 
-  // configure the rotor parameters
+  // process rotor parameters
-  Configure(rotor_element);
+  engine_power_est = Configure(rotor_element);
 
-  // configure the transmission parameters
+  // setup transmission if needed
   if (!ExternalRPM) {
-    Transmission = new FGTransmission(exec, num);
+
+    Transmission = new FGTransmission(exec, num, dt);
+
+    Transmission->SetThrusterMoment(PolarMoment);
+
+    // The MOI sensed behind the gear ( MOI_engine*sqr(GearRatio) ).
+    GearMoment = ConfigValueConv(rotor_element, "gearmoment", 0.1*PolarMoment, "SLUG*FT2");
+    GearMoment = Constrain(1e-6, GearMoment, 1e9);
+    Transmission->SetEngineMoment(GearMoment);
 
     Transmission->SetMaxBrakePower(MaxBrakePower);
 
-    if (rotor_element->FindElement("gearloss")) {
+    GearLoss = ConfigValueConv(rotor_element, "gearloss", 0.0025 * engine_power_est, "HP");
-      GearLoss = rotor_element->FindElementValueAsNumberConvertTo("gearloss","HP");
+    GearLoss = Constrain(0.0, GearLoss, 1e9);
-      GearLoss *= hptoftlbssec;
+    GearLoss *= hptoftlbssec;
-    }
-
-    if (GearLoss<1e-6) { // TODO, allow 0 ?
-      double ehp = 0.5 * BladeNum*BladeChord*Radius*Radius; // guess engine power
-      //cout << "# guessed engine power: " << ehp << endl;
-      GearLoss = 0.0025 * ehp * hptoftlbssec;
-    }
     Transmission->SetEngineFriction(GearLoss);
 
-    // The MOI sensed behind the gear ( MOI_engine*sqr(GearRatio) ).
-    if (rotor_element->FindElement("gearmoment")) {
-      GearMoment = rotor_element->FindElementValueAsNumberConvertTo("gearmoment","SLUG*FT2");
-    }
-
-    if (GearMoment<1e-2) { // TODO, need better check for lower limit
-      GearMoment = 0.1*PolarMoment;
-    }
-    Transmission->SetEngineMoment(GearMoment);
-
-    Transmission->SetThrusterMoment(PolarMoment);
   }
 
-  // shaft representation - a rather simple transform, 
+  // shaft representation - a rather simple transform,
   // but using a matrix is safer.
   TboToHsr.InitMatrix(   0.0, 0.0, 1.0,
                          0.0, 1.0, 0.0,
@@ -342,12 +223,11 @@ FGRotor::~FGRotor(){
   Debug(1);
 }
 
-
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-// 5in1: value-fetch-convert-default-return function 
+// 5in1: value-fetch-convert-default-return function
 
-double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val, 
+double FGRotor::ConfigValueConv( Element* el, const string& ename, double default_val,
                                   const string& unit, bool tell)
 {
 
@@ -388,20 +268,21 @@ double FGRotor::ConfigValue(Element* el, const string& ename, double default_val
 
 // 1. read configuration and try to fill holes, ymmv
 // 2. calculate derived parameters
-void FGRotor::Configure(Element* rotor_element)
+double FGRotor::Configure(Element* rotor_element)
 {
 
-  double estimate;
+  double estimate, engine_power_est=0.0;
   const bool yell   = true;
   const bool silent = false;
 
 
-  Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell); 
+  Radius = 0.5 * ConfigValueConv(rotor_element, "diameter", 42.0, "FT", yell);
   Radius = Constrain(1e-3, Radius, 1e9);
   
   BladeNum = (int) ConfigValue(rotor_element, "numblades", 3 , yell);
   
   GearRatio = ConfigValue(rotor_element, "gearratio", 1.0, yell);
+  GearRatio = Constrain(1e-9, GearRatio, 1e9);
 
   // make sure that v_tip (omega*r) is below 0.7mach ~ 750ft/s
   estimate = (750.0/Radius)/(2.0*M_PI) * 60.0;  // 7160/Radius
@@ -425,35 +306,28 @@ void FGRotor::Configure(Element* rotor_element)
 
   estimate = sqr(BladeChord) * sqr(Radius - HingeOffset) * 0.57;
   BladeFlappingMoment = ConfigValueConv(rotor_element, "flappingmoment", estimate, "SLUG*FT2");   
-  BladeFlappingMoment = Constrain(1.0e-6, BladeFlappingMoment, 1e9);
+  BladeFlappingMoment = Constrain(1e-9, BladeFlappingMoment, 1e9);
 
   // guess mass from moment of a thin stick, and multiply by the blades cg distance
   estimate = ( 3.0 * BladeFlappingMoment / sqr(Radius) ) * (0.45 * Radius) ;
   BladeMassMoment = ConfigValue(rotor_element, "massmoment", estimate); // unit is slug-ft
-  BladeMassMoment = Constrain(0.001, BladeMassMoment, 1e9);
+  BladeMassMoment = Constrain(1e-9, BladeMassMoment, 1e9);
 
   estimate = 1.1 * BladeFlappingMoment * BladeNum;
   PolarMoment = ConfigValueConv(rotor_element, "polarmoment", estimate, "SLUG*FT2");
-  PolarMoment = Constrain(1e-6, PolarMoment, 1e9);
+  PolarMoment = Constrain(1e-9, PolarMoment, 1e9);
 
   // "inflowlag" is treated further down.
 
   TipLossB = ConfigValue(rotor_element, "tiplossfactor", 1.0, silent);
 
-  estimate = 0.01 * PolarMoment ; // guesses for huey, bo105 20-30hp
+  // estimate engine power (bit of a pity, cause our caller already knows)
+  engine_power_est = 0.5 * BladeNum*BladeChord*Radius*Radius;
+
+  estimate = engine_power_est / 30.0;
   MaxBrakePower  = ConfigValueConv(rotor_element, "maxbrakepower", estimate, "HP");
   MaxBrakePower *= hptoftlbssec;
 
-  // ground effect
-  if (rotor_element->FindElement("cgroundeffect")) {
-    double cge,gee;
-    cge = rotor_element->FindElementValueAsNumber("cgroundeffect");
-    cge = Constrain(1e-9, cge, 1.0);
-    gee = 1.0 / ( 2.0*Radius * cge );
-    cerr << "# *** 'cgroundeffect' is defunct." << endl;
-    cerr << "# *** use 'groundeffectexp' with: " << gee << endl;
-  }
-
   GroundEffectExp = ConfigValue(rotor_element, "groundeffectexp", 0.0);
   GroundEffectShift = ConfigValueConv(rotor_element, "groundeffectshift", 0.0, "FT");
 
@@ -470,9 +344,9 @@ void FGRotor::Configure(Element* rotor_element)
   estimate = 16.0/(LockNumberByRho*rho * omega_tmp ); // 16/(gamma*Omega)
   // printf("# Est. InflowLag: %f\n", estimate);
   InflowLag = ConfigValue(rotor_element, "inflowlag", estimate, yell);
-  InflowLag = Constrain(1.0e-6, InflowLag, 2.0);
+  InflowLag = Constrain(1e-6, InflowLag, 2.0);
 
-  return;
+  return engine_power_est;
 } // Configure
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -507,10 +381,10 @@ FGColumnVector3 FGRotor::hub_vel_body2ca( const FGColumnVector3 &uvw,
 
 FGColumnVector3 FGRotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
 {
-  FGColumnVector3 av_s_fus, av_w_fus;    
+  FGColumnVector3 av_s_fus, av_w_fus;
 
   // for comparison:
-  // av_s_fus = BodyToShaft * pqr; /SH79/ 
+  // av_s_fus = BodyToShaft * pqr; /SH79/
   // BodyToShaft = TboToHsr * InvTransform
   av_s_fus = TboToHsr * InvTransform * pqr;
 
@@ -531,7 +405,7 @@ FGColumnVector3 FGRotor::fus_angvel_body2ca( const FGColumnVector3 &pqr)
 // reduction of inflow if the helicopter is close to the ground, yielding to
 // higher thrust, see /TA77/ eqn(10a).
 
-void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww, 
+void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww,
                                     double flow_scale)
 {
 
@@ -575,8 +449,8 @@ void FGRotor::calc_flow_and_thrust( double theta_0, double Uw, double Ww,
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-// The coning angle doesn't apply for teetering rotors, but calculating
+// Two blade teetering rotors are often 'preconed' to a fixed angle, but the 
-// doesn't hurt. /SH79/ eqn(29)
+// calculated value is pretty close to the real one. /SH79/ eqn(29)
 
 void FGRotor::calc_coning_angle(double theta_0)
 {
@@ -657,7 +531,7 @@ void FGRotor::calc_torque(double theta_0)
   // estimate blade drag
   double delta_dr = 0.009 + 0.3*sqr(6.0*C_T/(LiftCurveSlope*Solidity));
 
-  Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] * 
+  Torque = rho * BladeNum * BladeChord * delta_dr * sqr(Omega*Radius) * R[2] *
            (1.0+4.5*sqr(mu))/8.0
                      - (Thrust*lambda + H_drag*mu)*Radius;
 
@@ -666,6 +540,25 @@ void FGRotor::calc_torque(double theta_0)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+// Get the downwash angles with respect to the shaft axis.
+// Given a 'regular'  main rotor, the angles are zero when the downwash points
+// down, positive theta values mean that the downwash turns towards the nose,
+// and positive phi values mean it turns to the left side. (Note: only airspeed
+// is transformed, the rotational speed contribution is ignored.)
+
+void FGRotor::calc_downwash_angles()
+{
+  FGColumnVector3 v_shaft;
+  v_shaft = TboToHsr * InvTransform * in.AeroUVW;
+
+  theta_downwash = atan2( -v_shaft(eU), v_induced - v_shaft(eW)) + a1s;
+  phi_downwash   = atan2(  v_shaft(eV), v_induced - v_shaft(eW)) + b1s;
+
+  return;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 // transform rotor forces from control axes to shaft axes, and express
 // in body axes /SH79/ eqn(40,41)
 
@@ -674,7 +567,7 @@ FGColumnVector3 FGRotor::body_forces(double a_ic, double b_ic)
   FGColumnVector3 F_s(
         - H_drag*cos(beta_orient) - J_side*sin(beta_orient) + Thrust*b_ic,
         - H_drag*sin(beta_orient) + J_side*cos(beta_orient) + Thrust*a_ic,
-        - Thrust);    
+        - Thrust);
 
   return HsrToTbo * F_s;
 }
@@ -718,6 +611,7 @@ void FGRotor::CalcRotorState(void)
   // fetch needed values from environment
   rho = in.Density; // slugs/ft^3.
   double h_agl_ft = in.H_agl;
+
   // update InvTransform, the rotor orientation could have been altered
   InvTransform = Transform().Transposed();
 
@@ -736,16 +630,18 @@ void FGRotor::CalcRotorState(void)
   B_IC      = LongitudinalCtrl;
   theta_col = CollectiveCtrl;
 
-  // ground effect
+  // optional ground effect, a ge_factor of 1.0 gives no effect
+  // and 0.5 yields to maximal influence.
   if (GroundEffectExp > 1e-5) {
     if (h_agl_ft<0.0) h_agl_ft = 0.0; // clamp
     filtered_hagl = damp_hagl.execute(h_agl_ft) + GroundEffectShift;
     // actual/nominal factor avoids absurd scales at startup
-    ge_factor -= exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM);
+    ge_factor -= GroundEffectScaleNorm *
-    if (ge_factor<0.5) ge_factor=0.5; // clamp
+                 ( exp(-filtered_hagl*GroundEffectExp) * (RPM / NominalRPM) );
+    ge_factor = Constrain(0.5, ge_factor, 1.0);
   }
 
-  // all set, start calculations
+  // all set, start calculations ...
 
   vHub_ca  = hub_vel_body2ca(in.AeroUVW, in.AeroPQR, A_IC, B_IC);
 
@@ -761,12 +657,12 @@ void FGRotor::CalcRotorState(void)
 
   calc_torque(theta_col);
 
-  // Fixme: only valid for a 'decent' rotor
+  calc_downwash_angles();
-  theta_downwash = atan2( -in.AeroUVW(eU), v_induced - in.AeroUVW(eW));
-  phi_downwash   = atan2(  in.AeroUVW(eV), v_induced - in.AeroUVW(eW));
 
+  // ... and assign to inherited vFn and vMn members
+  //     (for processing see FGForce::GetBodyForces).
   vFn = body_forces(A_IC, B_IC);
-  vMn = Transform() * body_moments(A_IC, B_IC); 
+  vMn = Transform() * body_moments(A_IC, B_IC);
 
 }
 
@@ -778,8 +674,6 @@ double FGRotor::Calculate(double EnginePower)
   CalcRotorState();
 
   if (! ExternalRPM) {
-    Transmission->SetClutchCtrlNorm(ClutchCtrlNorm);
-
     // the RPM values are handled inside Transmission
     Transmission->Calculate(EnginePower, Torque, in.TotalDeltaT);
 
@@ -809,9 +703,6 @@ bool FGRotor::BindModel(void)
   property_name = base_property_name + "/engine-rpm";
   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetEngineRPM );
 
-  property_name = base_property_name + "/rotor-thrust-lbs"; // might be redundant - check!
-  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetThrust );
-
   property_name = base_property_name + "/a0-rad";
   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetA0 );
 
@@ -845,6 +736,10 @@ bool FGRotor::BindModel(void)
   property_name = base_property_name + "/phi-downwash-rad";
   PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetPhiDW );
 
+  property_name = base_property_name + "/groundeffect-scale-norm";
+  PropertyManager->Tie( property_name.c_str(), this, &FGRotor::GetGroundEffectScaleNorm,
+                                                     &FGRotor::SetGroundEffectScaleNorm );
+
   switch (ControlMap) {
     case eTailCtrl:
       property_name = base_property_name + "/antitorque-ctrl-rad";
@@ -891,7 +786,7 @@ bool FGRotor::BindModel(void)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-string FGRotor::GetThrusterLabels(int id, string delimeter)
+string FGRotor::GetThrusterLabels(int id, const string& delimeter)
 {
 
   ostringstream buf;
@@ -904,7 +799,7 @@ string FGRotor::GetThrusterLabels(int id, string delimeter)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-string FGRotor::GetThrusterValues(int id, string delimeter)
+string FGRotor::GetThrusterValues(int id, const string& delimeter)
 {
 
   ostringstream buf;
@@ -1003,5 +898,5 @@ void FGRotor::Debug(int from)
 }
 
 
-} // namespace JSBSim 
+} // namespace JSBSim
 

src/FDM/JSBSim/models/propulsion/FGRotor.h

@@ -28,6 +28,7 @@ HISTORY
 01/01/10  T.Kreitler test implementation
 01/10/11  T.Kreitler changed to single rotor model
 03/06/11  T.Kreitler added brake, clutch, and experimental free-wheeling-unit
+02/05/12  T.Kreitler brake, clutch, and FWU now in FGTransmission class
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 SENTRY
@@ -41,12 +42,13 @@ INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 #include "FGThruster.h"
+#include "FGTransmission.h"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_ROTOR "$Id: FGRotor.h,v 1.12 2011/10/15 21:30:28 jentron Exp $"
+#define ID_ROTOR "$Id: FGRotor.h,v 1.14 2012/03/18 15:48:36 jentron Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -89,7 +91,6 @@ CLASS DOCUMENTATION
   <groundeffectexp> {number} </groundeffectexp>
   <groundeffectshift unit="{LENGTH}"> {number} </groundeffectshift>
 
-  <freewheelthresh> {number} </freewheelthresh>
 </rotor>
 
 //  LENGTH means any of the supported units, same for ANGLE and MOMENT.
@@ -135,28 +136,29 @@ CLASS DOCUMENTATION
     Experimental properties
 
     \<groundeffectexp>    - Exponent for ground effect approximation. Values usually range from 0.04
-                            for large rotors to 0.1 for smaller ones. As a rule of thumb the effect 
+                            for large rotors to 0.1 for smaller ones. As a rule of thumb the effect
                             vanishes at a height 2-3 times the rotor diameter.
                               formula used: exp ( - groundeffectexp * (height+groundeffectshift) )
                             Omitting or setting to 0.0 disables the effect calculation.
-    \<groundeffectshift>  - Further adjustment of ground effect, approx. hub height or slightly above. 
+    \<groundeffectshift>  - Further adjustment of ground effect, approx. hub height or slightly above
+                            (This lessens the influence of the ground effect).
 
 </pre>
 
-<h3>Notes:</h3>  
+<h3>Notes:</h3>
 
   <h4>- Controls -</h4>
 
     The behavior of the rotor is controlled/influenced by following inputs.<ul>
       <li> The power provided by the engine. This is handled by the regular engine controls.</li>
-      <li> The collective control input. This is read from the <tt>fdm</tt> property 
+      <li> The collective control input. This is read from the <tt>fdm</tt> property
            <tt>propulsion/engine[x]/collective-ctrl-rad</tt>. See below for tail rotor</li>
       <li> The lateral cyclic input. Read from
            <tt>propulsion/engine[x]/lateral-ctrl-rad</tt>.</li>
       <li> The longitudinal cyclic input. Read from 
            <tt>propulsion/engine[x]/longitudinal-ctrl-rad</tt>.</li>
-      <li> The tail collective (aka antitorque, aka pedal) control input. Read from
+      <li> The tail rotor collective (aka antitorque, aka pedal) control input. Read from
-           <tt>propulsion/engine[x]/antitorque-ctrl-rad</tt> or 
+           <tt>propulsion/engine[x]/antitorque-ctrl-rad</tt> or
            <tt>propulsion/engine[x]/tail-collective-ctrl-rad</tt>.</li>
 
     </ul>
@@ -167,7 +169,7 @@ CLASS DOCUMENTATION
     is linked to to the main (=first, =0) rotor, and specifing
     <tt>\<controlmap\> TAIL \</controlmap\></tt> tells this rotor to read the
     collective input from <tt>propulsion/engine[1]/antitorque-ctrl-rad</tt>
-    (The TAIL-map ignores lateral and longitudinal input). The rotor needs to be 
+    (The TAIL-map ignores lateral and longitudinal input). The rotor needs to be
     attached to a dummy engine, e.g. an 1HP electrical engine.
     A tandem rotor is setup analogous. 
 
@@ -198,6 +200,12 @@ CLASS DOCUMENTATION
 
     </ul>
 
+  <h4>- Scaling the ground effect -</h4>
+
+    The property <tt>propulsion/engine[x]/groundeffect-scale-norm</tt> allows fdm based
+    scaling of the ground effect influence. For instance the effect vanishes at speeds
+    above approx. 50kts, or one likes to land on a 'perforated' helipad.
+
   <h4>- Development hints -</h4>
 
     Setting <tt>\<ExternalRPM> -1 \</ExternalRPM></tt> the rotor's RPM is controlled  by
@@ -205,24 +213,24 @@ CLASS DOCUMENTATION
     when developing a FDM.
   
 
-<h3>References:</h3>  
+<h3>References:</h3>
 
     <dl>    
     <dt>/SH79/</dt><dd>Shaugnessy, J. D., Deaux, Thomas N., and Yenni, Kenneth R.,
-              "Development and Validation of a Piloted Simulation of a 
+              "Development and Validation of a Piloted Simulation of a
               Helicopter and External Sling Load",  NASA TP-1285, 1979.</dd>
     <dt>/BA41/</dt><dd>Bailey,F.J.,Jr., "A Simplified Theoretical Method of Determining
               the Characteristics of a Lifting Rotor in Forward Flight", NACA Rep.716, 1941.</dd>
     <dt>/AM50/</dt><dd>Amer, Kenneth B.,"Theory of Helicopter Damping in Pitch or Roll and a
               Comparison With Flight Measurements", NACA TN-2136, 1950.</dd>
     <dt>/TA77/</dt><dd>Talbot, Peter D., Corliss, Lloyd D., "A Mathematical Force and Moment
-              Model of a UH-1H Helicopter for Flight Dynamics Simulations", NASA TM-73,254, 1977.</dd> 
+              Model of a UH-1H Helicopter for Flight Dynamics Simulations", NASA TM-73,254, 1977.</dd>
     <dt>/GE49/</dt><dd>Gessow, Alfred, Amer, Kenneth B. "An Introduction to the Physical 
-              Aspects of Helicopter Stability", NACA TN-1982, 1949.</dd>  
+              Aspects of Helicopter Stability", NACA TN-1982, 1949.</dd>
     </dl>
 
     @author Thomas Kreitler
-    @version $Id: FGRotor.h,v 1.12 2011/10/15 21:30:28 jentron Exp $
+    @version $Id: FGRotor.h,v 1.14 2012/03/18 15:48:36 jentron Exp $
   */
 
 
@@ -231,57 +239,6 @@ CLASS DOCUMENTATION
 CLASS DECLARATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-class FGTransmission :  public FGJSBBase {
-
-public:
-  FGTransmission(FGFDMExec *exec, int num);
-  ~FGTransmission();
-
-  void Calculate(double EnginePower, double ThrusterTorque, double dt);
-
-  void   SetMaxBrakePower(double x) {MaxBrakePower=x;}
-  double GetMaxBrakePower() const {return MaxBrakePower;}
-  void   SetEngineFriction(double x) {EngineFriction=x;}
-  double GetEngineFriction() const {return EngineFriction;}
-  void   SetEngineMoment(double x) {EngineMoment=x;}
-  double GetEngineMoment() const {return EngineMoment;}
-  void   SetThrusterMoment(double x) {ThrusterMoment=x;}
-  double GetThrusterMoment() const {return ThrusterMoment;}
-
-  double GetFreeWheelTransmission() const {return FreeWheelTransmission;}
-  double GetEngineRPM() {return EngineRPM;}
-  double GetThrusterRPM() {return ThrusterRPM;}
-
-  double GetBrakeCtrl() const {return BrakeCtrlNorm;}
-  void   SetBrakeCtrl(double x) {BrakeCtrlNorm=x;}
-  void   SetClutchCtrlNorm(double x) {ClutchCtrlNorm=x;}
-
-private:
-  bool BindModel(int num);
-  // void Debug(int from);
-
-  inline double omega_to_rpm(double w) {return w * 9.54929658551372014613302580235;} // omega/(2.0*PI) * 60.0
-  inline double rpm_to_omega(double r) {return r * .104719755119659774615421446109;} // (rpm/60.0)*2.0*PI
-
-  Filter FreeWheelLag;
-  double FreeWheelTransmission; // state, 0: free, 1:locked
-
-  double ThrusterMoment;
-  double EngineMoment;   // estimated MOI of gear and engine, influences acceleration
-  double EngineFriction; // estimated friction in gear and possibly engine
-
-  double ClutchCtrlNorm; // also in FGThruster.h
-  double BrakeCtrlNorm;
-  double MaxBrakePower;
-
-  double EngineRPM;
-  double ThrusterRPM;
-  FGPropertyManager* PropertyManager;
-
-};
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 class FGRotor :  public FGThruster {
 
   enum eCtrlMapping {eMainCtrl=0, eTailCtrl, eTandemCtrl};
@@ -297,10 +254,10 @@ public:
   /// Destructor for FGRotor
   ~FGRotor();
 
-  /** Returns the power required by the rotor. */
+  /// Returns the power required by the rotor.
   double GetPowerRequired(void)const { return PowerRequired; }
 
-  /** Returns the scalar thrust of the rotor, and adjusts the RPM value. */
+  /// Returns the scalar thrust of the rotor, and adjusts the RPM value.
   double Calculate(double EnginePower);
 
 
@@ -336,11 +293,16 @@ public:
   /// Retrieves the torque
   double GetTorque(void) const { return Torque; }
   
-  /// Downwash angle - currently only valid for a rotor that spins horizontally
+  /// Downwash angle - positive values point forward (given a horizontal spinning rotor)
   double GetThetaDW(void) const { return theta_downwash; }
-  /// Downwash angle - currently only valid for a rotor that spins horizontally
+  /// Downwash angle - positive values point leftward (given a horizontal spinning rotor)
   double GetPhiDW(void) const { return phi_downwash; }
 
+  /// Retrieves the ground effect scaling factor.
+  double GetGroundEffectScaleNorm(void) const { return GroundEffectScaleNorm; }
+  /// Sets the ground effect scaling factor.
+  void   SetGroundEffectScaleNorm(double g) { GroundEffectScaleNorm = g; }
+
   /// Retrieves the collective control input in radians.
   double GetCollectiveCtrl(void) const { return CollectiveCtrl; }
   /// Retrieves the lateral control input in radians.
@@ -356,8 +318,8 @@ public:
   void SetLongitudinalCtrl(double c) { LongitudinalCtrl = c; }
 
   // Stubs. Only main rotor RPM is returned
-  string GetThrusterLabels(int id, string delimeter);
+  string GetThrusterLabels(int id, const string& delimeter);
-  string GetThrusterValues(int id, string delimeter);
+  string GetThrusterValues(int id, const string& delimeter);
 
 private:
 
@@ -368,7 +330,7 @@ private:
   double ConfigValue( Element* e, const string& ename, double default_val=0.0,
                                   bool tell=false);
 
-  void Configure(Element* rotor_element);
+  double Configure(Element* rotor_element);
 
   void CalcRotorState(void);
 
@@ -378,6 +340,7 @@ private:
   void calc_flapping_angles(double theta_0, const FGColumnVector3 &pqr_fus_w);
   void calc_drag_and_side_forces(double theta_0);
   void calc_torque(double theta_0);
+  void calc_downwash_angles();
 
   // transformations
   FGColumnVector3 hub_vel_body2ca( const FGColumnVector3 &uvw, const FGColumnVector3 &pqr, 
@@ -409,6 +372,7 @@ private:
   FGPropertyManager* ExtRPMsource;
   double SourceGearRatio;
 
+  // 'real' rotor parameters
   double BladeChord;
   double LiftCurveSlope;
   double BladeTwist;
@@ -419,8 +383,10 @@ private:
   double InflowLag;
   double TipLossB;
 
+  // groundeffect
   double GroundEffectExp;
   double GroundEffectShift;
+  double GroundEffectScaleNorm;
 
   // derived parameters
   double LockNumberByRho;
@@ -449,7 +415,7 @@ private:
   double lambda;     // inflow ratio
   double mu;         // tip-speed ratio 
   double nu;         // induced inflow ratio
-  double v_induced;  // induced velocity, always positive [ft/s]
+  double v_induced;  // induced velocity, usually positive [ft/s]
 
   double theta_downwash;
   double phi_downwash;

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

@@ -45,7 +45,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id: FGThruster.cpp,v 1.14 2011/03/10 01:35:25 dpculp Exp $";
+static const char *IdSrc = "$Id: FGThruster.cpp,v 1.16 2012/03/18 15:48:36 jentron Exp $";
 static const char *IdHdr = ID_THRUSTER;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -66,7 +66,6 @@ FGThruster::FGThruster(FGFDMExec *FDMExec, Element *el, int num ): FGForce(FDMEx
 
   GearRatio = 1.0;
   ReverserAngle = 0.0;
-  ClutchCtrlNorm = 1.0;
   EngineNum = num;
   PropertyManager = FDMExec->GetPropertyManager();
 
@@ -99,13 +98,6 @@ FGThruster::FGThruster(FGFDMExec *FDMExec, Element *el, int num ): FGForce(FDMEx
                                                           &FGThruster::SetReverserAngle);
   }
 
-  if (el->GetName() == "rotor") // At this time only a rotor can have a clutch.
-  {
-    property_name = base_property_name + "/clutch-ctrl-norm";
-    PropertyManager->Tie( property_name.c_str(), (FGThruster *)this, &FGThruster::GetClutchCtrl,
-                                                          &FGThruster::SetClutchCtrl);
-  }
-
   Debug(0);
 }
 
@@ -118,7 +110,7 @@ FGThruster::~FGThruster()
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-string FGThruster::GetThrusterLabels(int id, string delimeter)
+string FGThruster::GetThrusterLabels(int id, const string& delimeter)
 {
   std::ostringstream buf;
 
@@ -129,7 +121,7 @@ string FGThruster::GetThrusterLabels(int id, string delimeter)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-string FGThruster::GetThrusterValues(int id, string delimeter)
+string FGThruster::GetThrusterValues(int id, const string& delimeter)
 {
   std::ostringstream buf;
 

src/FDM/JSBSim/models/propulsion/FGThruster.h

@@ -46,7 +46,7 @@ INCLUDES
 DEFINITIONS
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-#define ID_THRUSTER "$Id: FGThruster.h,v 1.18 2011/09/24 14:26:46 jentron Exp $"
+#define ID_THRUSTER "$Id: FGThruster.h,v 1.20 2012/03/18 15:48:36 jentron Exp $"
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FORWARD DECLARATIONS
@@ -74,7 +74,7 @@ CLASS DOCUMENTATION
     1.57 (pi/2) results in no thrust at all.
  
     @author Jon Berndt
-    @version $Id: FGThruster.h,v 1.18 2011/09/24 14:26:46 jentron Exp $
+    @version $Id: FGThruster.h,v 1.20 2012/03/18 15:48:36 jentron Exp $
     */
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -106,13 +106,11 @@ public:
   string GetName(void) {return Name;}
   void SetReverserAngle(double angle) {ReverserAngle = angle;}
   double GetReverserAngle(void) const {return ReverserAngle;}
-  double GetClutchCtrl(void) const { return ClutchCtrlNorm; }
-  void SetClutchCtrl(double c) { ClutchCtrlNorm = c; }
   virtual double GetRPM(void) const { return 0.0; };
   virtual double GetEngineRPM(void) const { return 0.0; };
   double GetGearRatio(void) {return GearRatio; }
-  virtual string GetThrusterLabels(int id, string delimeter);
+  virtual string GetThrusterLabels(int id, const string& delimeter);
-  virtual string GetThrusterValues(int id, string delimeter);
+  virtual string GetThrusterValues(int id, const string& delimeter);
 
   struct Inputs {
     double TotalDeltaT;
@@ -137,7 +135,6 @@ protected:
   double GearRatio;
   double ThrustCoeff;
   double ReverserAngle;
-  double ClutchCtrlNorm;  
   int EngineNum;
   FGPropertyManager* PropertyManager;
   virtual void Debug(int from);

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

@@ -0,0 +1,206 @@
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ JSBSim
+ Author:       Jon S. Berndt
+ Date started: 08/24/00
+ ------------- Copyright (C) 2000  Jon S. Berndt (jon@jsbsim.org) -------------
+
+ Module:       FGTransmission.cpp
+ Author:       T.Kreitler
+ Date started: 02/05/12
+
+ ------------- Copyright (C) 2012  T. Kreitler (t.kreitler@web.de) -------------
+
+ This program is free software; you can redistribute it and/or modify it under
+ 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 Lesser General Public License for more
+ details.
+
+ 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 Lesser General Public License can also be found on
+ the world wide web at http://www.gnu.org.
+
+FUNCTIONAL DESCRIPTION
+--------------------------------------------------------------------------------
+
+HISTORY
+--------------------------------------------------------------------------------
+02/05/12  T.Kreitler  Created
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+INCLUDES
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+
+#include "FGTransmission.h"
+
+using std::cout;
+using std::endl;
+
+namespace JSBSim {
+
+static const char *IdSrc = "$Id: FGTransmission.cpp,v 1.1 2012/02/25 14:37:02 jentron Exp $";
+static const char *IdHdr = ID_TRANSMISSION;
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CLASS IMPLEMENTATION
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+FGTransmission::FGTransmission(FGFDMExec *exec, int num, double dt) :
+  FreeWheelTransmission(1.0),
+  ThrusterMoment(1.0), EngineMoment(1.0), EngineFriction(0.0),
+  ClutchCtrlNorm(1.0), BrakeCtrlNorm(0.0), MaxBrakePower(0.0),
+  EngineRPM(0.0), ThrusterRPM(0.0)
+{
+  PropertyManager = exec->GetPropertyManager();
+  FreeWheelLag = Filter(200.0,dt); // avoid too abrupt changes in transmission
+  BindModel(num);
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+FGTransmission::~FGTransmission(){
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+// basically P = Q*w and Q_Engine + (-Q_Rotor) = J * dw/dt, J = Moment
+//
+void FGTransmission::Calculate(double EnginePower, double ThrusterTorque, double dt) {
+
+  double coupling = 1.0, coupling_sq = 1.0;
+  double fw_mult = 1.0;
+
+  double d_omega = 0.0, engine_d_omega = 0.0, thruster_d_omega = 0.0; // relative changes
+
+  double engine_omega = rpm_to_omega(EngineRPM);
+  double safe_engine_omega = engine_omega < 1e-1 ? 1e-1 : engine_omega;
+  double engine_torque = EnginePower / safe_engine_omega;
+
+  double thruster_omega = rpm_to_omega(ThrusterRPM);
+  double safe_thruster_omega = thruster_omega < 1e-1 ? 1e-1 : thruster_omega;
+
+  engine_torque  -= EngineFriction / safe_engine_omega;
+  ThrusterTorque += Constrain(0.0, BrakeCtrlNorm, 1.0) * MaxBrakePower / safe_thruster_omega;
+
+  // would the FWU release ?
+  engine_d_omega = engine_torque/EngineMoment * dt;
+  thruster_d_omega =  - ThrusterTorque/ThrusterMoment * dt;
+
+  if ( thruster_omega+thruster_d_omega > engine_omega+engine_d_omega ) {
+    // don't drive the engine
+    FreeWheelTransmission = 0.0;
+  } else {
+    FreeWheelTransmission = 1.0;
+  }
+
+  fw_mult = FreeWheelLag.execute(FreeWheelTransmission);
+  coupling = fw_mult * Constrain(0.0, ClutchCtrlNorm, 1.0);
+
+  if (coupling < 0.999999) { // are the separate calculations needed ?
+    // assume linear transfer 
+    engine_d_omega   =
+       (engine_torque - ThrusterTorque*coupling)/(ThrusterMoment*coupling + EngineMoment) * dt;
+    thruster_d_omega = 
+       (engine_torque*coupling - ThrusterTorque)/(ThrusterMoment + EngineMoment*coupling) * dt;
+
+    EngineRPM += omega_to_rpm(engine_d_omega);
+    ThrusterRPM += omega_to_rpm(thruster_d_omega);
+
+    // simulate transit to static friction
+    coupling_sq = coupling*coupling;
+    EngineRPM   = (1.0-coupling_sq) * EngineRPM    + coupling_sq * 0.02 * (49.0*EngineRPM + ThrusterRPM);
+    ThrusterRPM = (1.0-coupling_sq) * ThrusterRPM  + coupling_sq * 0.02 * (EngineRPM + 49.0*ThrusterRPM);
+
+    // enforce equal rpm
+    if ( fabs(EngineRPM-ThrusterRPM) < 1e-3 ) {
+      EngineRPM = ThrusterRPM = 0.5 * (EngineRPM+ThrusterRPM);
+    }
+  } else {
+    d_omega = (engine_torque - ThrusterTorque)/(ThrusterMoment + EngineMoment) * dt;
+    EngineRPM = ThrusterRPM += omega_to_rpm(d_omega);
+  }
+
+  // nothing will turn backward
+  if (EngineRPM < 0.0 ) EngineRPM = 0.0;
+  if (ThrusterRPM < 0.0 ) ThrusterRPM = 0.0;
+
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+bool FGTransmission::BindModel(int num)
+{
+  string property_name, base_property_name;
+  base_property_name = CreateIndexedPropertyName("propulsion/engine", num);
+
+  property_name = base_property_name + "/brake-ctrl-norm";
+  PropertyManager->Tie( property_name.c_str(), this, 
+      &FGTransmission::GetBrakeCtrlNorm, &FGTransmission::SetBrakeCtrlNorm);
+
+  property_name = base_property_name + "/clutch-ctrl-norm";
+  PropertyManager->Tie( property_name.c_str(), this, 
+      &FGTransmission::GetClutchCtrlNorm, &FGTransmission::SetClutchCtrlNorm);
+
+  property_name = base_property_name + "/free-wheel-transmission";
+  PropertyManager->Tie( property_name.c_str(), this, 
+      &FGTransmission::GetFreeWheelTransmission);
+
+  return true;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//    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 FGTransmission::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: FGTransmission" << endl;
+    if (from == 1) cout << "Destroyed:    FGTransmission" << 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 & 64) {
+    if (from == 0) { // Constructor
+      cout << IdSrc << endl;
+      cout << IdHdr << endl;
+    }
+  }
+}
+
+}
+

src/FDM/JSBSim/models/propulsion/FGTransmission.h

@@ -0,0 +1,179 @@
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ JSBSim
+ Author:       Jon S. Berndt
+ Date started: 08/24/00
+ ------------- Copyright (C) 2000  Jon S. Berndt (jon@jsbsim.org) -------------
+
+ Header:       FGTransmission.h
+ Author:       T.Kreitler
+ Date started: 02/05/12
+
+ ------------- Copyright (C) 2012  T. Kreitler (t.kreitler@web.de) -------------
+
+ This program is free software; you can redistribute it and/or modify it under
+ 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 Lesser General Public License for more
+ details.
+
+ 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 Lesser General Public License can also be found on
+ the world wide web at http://www.gnu.org.
+
+HISTORY
+--------------------------------------------------------------------------------
+02/05/12  T.Kreitler  Created
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+SENTRY
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+#ifndef FGTRANSMISSION_H
+#define FGTRANSMISSION_H
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+INCLUDES
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+#include "FGJSBBase.h"
+#include "FGFDMExec.h"
+#include "input_output/FGPropertyManager.h"
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+DEFINITIONS
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+#define ID_TRANSMISSION "$Id: FGTransmission.h,v 1.1 2012/02/25 14:37:02 jentron Exp $"
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+FORWARD DECLARATIONS
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+namespace JSBSim {
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CLASS DOCUMENTATION
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+/** Utility class that handles power transmission in conjunction with FGRotor.
+
+  This class provides brake, clutch and free-wheel-unit (FWU) functionality
+  for the rotor model. Also it is responsible for the RPM calculations.
+
+  When the engine is off the brake could be used to slow/hold down a spinning
+  rotor. The maximum brake power is defined in the rotors' config file.
+  (Right now there is no checking if the input is in the [0..1] range.)
+
+  The clutch operation is based on a heuristic approach. In the intermediate
+  state the transfer is proportional to the clutch position. But equal RPM
+  values are enforced on the thruster and rotor sides when approaching the
+  closed state.
+
+  The FWU inhibits that the rotor is driving the engine. To do so, the code
+  just predicts the upcoming FWU state based on current torque conditions.
+
+  Some engines won't work properly when the clutch is open. To keep them
+  controllable some load must be provided on the engine side (EngineFriction,
+  aka gear-loss). See the notes under 'Engine issues' in FGRotor.
+
+<h3>Property tree</h3>
+
+  The following properties are created (with x = your thruster number):
+  <pre>
+    propulsion/engine[x]/brake-ctrl-norm
+    propulsion/engine[x]/free-wheel-transmission
+    propulsion/engine[x]/clutch-ctrl-norm
+  </pre>
+
+<h3>Notes</h3>
+
+  <ul>
+    <li> EngineFriction is assumed constant, so better orientate at low RPM
+         values, because piston and turboprop engines don't 'like' high
+         load at startup.</li>
+    <li> The model doesn't support backward operation.</li>
+    <li> And even worse, the torque calculations silently assume a minimal
+         RPM value of approx. 1.</li>
+  </ul>
+
+*/
+
+
+/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CLASS DECLARATION
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+
+class FGTransmission : public FGJSBBase {
+
+public:
+  /** Constructor for FGTransmission.
+      @param exec a pointer to the main executive object
+      @param num the number of the thruster that uses this object
+      @param dt simulation delta T */
+  FGTransmission(FGFDMExec *exec, int num, double dt);
+
+  /// Destructor for FGTransmission
+  ~FGTransmission();
+
+  void Calculate(double EnginePower, double ThrusterTorque, double dt);
+
+  void   SetMaxBrakePower(double x) {MaxBrakePower=x;}
+  double GetMaxBrakePower() const {return MaxBrakePower;}
+  void   SetEngineFriction(double x) {EngineFriction=x;}
+  double GetEngineFriction() const {return EngineFriction;}
+  void   SetEngineMoment(double x) {EngineMoment=x;}
+  double GetEngineMoment() const {return EngineMoment;}
+  void   SetThrusterMoment(double x) {ThrusterMoment=x;}
+  double GetThrusterMoment() const {return ThrusterMoment;}
+
+  double GetFreeWheelTransmission() const {return FreeWheelTransmission;}
+  void   SetEngineRPM(double x) {EngineRPM=x;}
+  double GetEngineRPM() {return EngineRPM;}
+  void   SetThrusterRPM(double x) {ThrusterRPM=x;}
+  double GetThrusterRPM() {return ThrusterRPM;}
+
+  double GetBrakeCtrlNorm() const {return BrakeCtrlNorm;}
+  void   SetBrakeCtrlNorm(double x) {BrakeCtrlNorm=x;}
+  double GetClutchCtrlNorm() const {return ClutchCtrlNorm;}
+  void   SetClutchCtrlNorm(double x) {ClutchCtrlNorm=x;}
+
+private:
+  bool BindModel(int num);
+  void Debug(int from);
+
+  inline double omega_to_rpm(double w) {
+    return w * 9.54929658551372014613302580235; // omega/(2.0*PI) * 60.0
+  }
+  inline double rpm_to_omega(double r) {
+    return r * 0.104719755119659774615421446109; // (rpm/60.0)*2.0*PI
+  }
+
+  Filter FreeWheelLag;
+  double FreeWheelTransmission; // state, 0: free, 1:locked
+
+  double ThrusterMoment;
+  double EngineMoment;   // estimated MOI of gear and engine, influences acceleration
+  double EngineFriction; // estimated friction in gear and possibly engine
+
+  double ClutchCtrlNorm;
+  double BrakeCtrlNorm;
+  double MaxBrakePower;
+
+  double EngineRPM;
+  double ThrusterRPM;
+  FGPropertyManager* PropertyManager;
+
+};
+
+}
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#endif
+