]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/JSBSim/models/FGLGear.cpp
Merge branch 'next' of http://git.gitorious.org/fg/flightgear into next
[flightgear.git] / src / FDM / JSBSim / models / FGLGear.cpp
index e8df4370ad99c828950097f6721c05ae415d02cb..4a4c9946f381e4ebd5aadc63f9ee2a10fdcc3441 100644 (file)
@@ -8,7 +8,7 @@
  Purpose:      Encapsulates the landing gear elements
  Called by:    FGAircraft
 
- ------------- Copyright (C) 1999  Jon S. Berndt (jsb@hal-pc.org) -------------
+ ------------- Copyright (C) 1999  Jon S. Berndt (jon@jsbsim.org) -------------
 
  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
@@ -41,6 +41,16 @@ INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 #include "FGLGear.h"
+#include "FGGroundReactions.h"
+#include "FGFCS.h"
+#include "FGAuxiliary.h"
+#include "FGAtmosphere.h"
+#include "FGMassBalance.h"
+#include "math/FGTable.h"
+#include <cstdlib>
+#include <cstring>
+
+using namespace std;
 
 namespace JSBSim {
 
@@ -52,16 +62,23 @@ DEFINITIONS
 GLOBAL DATA
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
-static const char *IdSrc = "$Id$";
+static const char *IdSrc = "$Id: FGLGear.cpp,v 1.80 2011/01/24 13:01:56 jberndt 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);
+
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CLASS IMPLEMENTATION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
+  FGForce(fdmex),
   GearNumber(number),
-  Exec(fdmex)
+  SteerAngle(0.0),
+  Castered(false),
+  StaticFriction(false)
 {
   Element *force_table=0;
   Element *dampCoeff=0;
@@ -81,7 +98,8 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
   } else if (sContactType == "STRUCTURE") {
     eContactType = ctSTRUCTURE;
   } else {
-    eContactType = ctUNKNOWN;
+    // Unknown contact point types will be treated as STRUCTURE.
+    eContactType = ctSTRUCTURE;
   }
 
   if (el->FindElement("spring_coeff"))
@@ -125,7 +143,7 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
   while (force_table) {
     force_type = force_table->GetAttributeValue("type");
     if (force_type == "CORNERING_COEFF") {
-      ForceY_Table = new FGTable(Exec->GetPropertyManager(), force_table);
+      ForceY_Table = new FGTable(fdmex->GetPropertyManager(), force_table);
     } else {
       cerr << "Undefined force table for " << name << " contact point" << endl;
     }
@@ -139,8 +157,37 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
   else sSteerType = "STEERABLE";
 
   Element* element = el->FindElement("location");
-  if (element) vXYZ = element->FindElementTripletConvertTo("IN");
+  if (element) vXYZn = element->FindElementTripletConvertTo("IN");
   else {cerr << "No location given for contact " << name << endl; exit(-1);}
+  SetTransformType(FGForce::tCustom);
+
+  element = el->FindElement("orientation");
+  if (element && (eContactType == ctBOGEY)) {
+    vGearOrient = element->FindElementTripletConvertTo("RAD");
+
+    double cp,sp,cr,sr,cy,sy;
+
+    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.;
+    mTGear(2,2) = 1.;
+    mTGear(3,3) = 1.;
+  }
 
   if      (sBrakeGroup == "LEFT"  ) eBrakeGrp = bgLeft;
   else if (sBrakeGroup == "RIGHT" ) eBrakeGrp = bgRight;
@@ -157,7 +204,7 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
 
   if      (sSteerType == "STEERABLE") eSteerType = stSteer;
   else if (sSteerType == "FIXED"    ) eSteerType = stFixed;
-  else if (sSteerType == "CASTERED" ) eSteerType = stCaster;
+  else if (sSteerType == "CASTERED" ) {eSteerType = stCaster; Castered = true;}
   else if (sSteerType.empty()       ) {eSteerType = stFixed;
                                        sSteerType = "FIXED (defaulted)";}
   else {
@@ -165,44 +212,11 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
          << sSteerType << " is undefined." << endl;
   }
 
-  RFRV = 0.7;  // Rolling force relaxation velocity, default value
-  SFRV = 0.7;  // Side force relaxation velocity, default value
-
-  Element* relax_vel = el->FindElement("relaxation_velocity");
-  if (relax_vel) {
-    if (relax_vel->FindElement("rolling")) {
-      RFRV = relax_vel->FindElementValueAsNumberConvertTo("rolling", "FT/SEC");
-    }
-    if (relax_vel->FindElement("side")) {
-      SFRV = relax_vel->FindElementValueAsNumberConvertTo("side", "FT/SEC");
-    }
-  }
-
-  State = Exec->GetState();
-  LongForceLagFilterCoeff = 1/State->Getdt(); // default longitudinal force filter coefficient
-  LatForceLagFilterCoeff  = 1/State->Getdt(); // default lateral force filter coefficient
-
-  Element* force_lag_filter_elem = el->FindElement("force_lag_filter");
-  if (force_lag_filter_elem) {
-    if (force_lag_filter_elem->FindElement("rolling")) {
-      LongForceLagFilterCoeff = force_lag_filter_elem->FindElementValueAsNumber("rolling");
-    }
-    if (force_lag_filter_elem->FindElement("side")) {
-      LatForceLagFilterCoeff = force_lag_filter_elem->FindElementValueAsNumber("side");
-    }
-  }
-
-  LongForceFilter = Filter(LongForceLagFilterCoeff, State->Getdt());
-  LatForceFilter = Filter(LatForceLagFilterCoeff, State->Getdt());
-
-  WheelSlipLagFilterCoeff = 1/State->Getdt();
-
-  Element *wheel_slip_angle_lag_elem = el->FindElement("wheel_slip_filter");
-  if (wheel_slip_angle_lag_elem) {
-    WheelSlipLagFilterCoeff = wheel_slip_angle_lag_elem->GetDataAsNumber();
-  }
-  
-  WheelSlipFilter = Filter(WheelSlipLagFilterCoeff, State->Getdt());
+  Auxiliary       = fdmex->GetAuxiliary();
+  Propagate       = fdmex->GetPropagate();
+  FCS             = fdmex->GetFCS();
+  MassBalance     = fdmex->GetMassBalance();
+  GroundReactions = fdmex->GetGroundReactions();
 
   GearUp = false;
   GearDown = true;
@@ -213,13 +227,6 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
 // Add some AI here to determine if gear is located properly according to its
 // brake group type ??
 
-  State       = Exec->GetState();
-  Aircraft    = Exec->GetAircraft();
-  Propagate   = Exec->GetPropagate();
-  Auxiliary   = Exec->GetAuxiliary();
-  FCS         = Exec->GetFCS();
-  MassBalance = Exec->GetMassBalance();
-
   WOW = lastWOW = false;
   ReportEnable = true;
   FirstContact = false;
@@ -229,11 +236,9 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
   MaximumStrutForce = MaximumStrutTravel = 0.0;
   SinkRate = GroundSpeed = 0.0;
 
-  vWhlBodyVec = MassBalance->StructuralToBody(vXYZ);
-
+  vWhlBodyVec = MassBalance->StructuralToBody(vXYZn);
   vLocalGear = Propagate->GetTb2l() * vWhlBodyVec;
-
-  vLocalWhlVel.InitMatrix();
+  vWhlVelVec.InitMatrix();
 
   compressLength  = 0.0;
   compressSpeed   = 0.0;
@@ -250,6 +255,9 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
   Peak = staticFCoeff;
   Curvature = 1.03;
 
+  // Initialize Lagrange multipliers
+  memset(LMultiplier, 0, sizeof(LMultiplier));
+
   Debug(0);
 }
 
@@ -263,36 +271,40 @@ FGLGear::~FGLGear()
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-FGColumnVector3& FGLGear::Force(void)
+FGColumnVector3& FGLGear::GetBodyForces(void)
 {
-  double t = Exec->GetState()->Getsim_time();
-  dT = State->Getdt()*Exec->GetGroundReactions()->GetRate();
+  double t = fdmex->GetSimTime();
+  dT = fdmex->GetDeltaT()*GroundReactions->GetRate();
 
-  vForce.InitMatrix();
-  vLocalForce.InitMatrix();
-  vMoment.InitMatrix();
+  vFn.InitMatrix();
 
   if (isRetractable) ComputeRetractionState();
 
   if (GearDown) {
+    FGColumnVector3 angularVel;
 
-    vWhlBodyVec = MassBalance->StructuralToBody(vXYZ); // Get wheel in body frame
+    vWhlBodyVec = MassBalance->StructuralToBody(vXYZn); // Get wheel in body frame
     vLocalGear = Propagate->GetTb2l() * vWhlBodyVec; // Get local frame wheel location
 
     gearLoc = Propagate->GetLocation().LocalToLocation(vLocalGear);
-    // Compute the height of the theoritical location of the wheel (if struct was not compressed) with
-    // respect to the ground level
-    double height = Exec->GetGroundCallback()->GetAGLevel(t, gearLoc, contact, normal, cvel);
-    vGroundNormal = -1. * Propagate->GetTec2b() * normal;
+    // Compute the height of the theoretical location of the wheel (if strut is
+    // not compressed) with respect to the ground level
+    double height = fdmex->GetGroundCallback()->GetAGLevel(t, gearLoc, contact, normal, cvel, angularVel);
+    vGroundNormal = Propagate->GetTec2b() * normal;
+
+    // The height returned above is the AGL and is expressed in the Z direction
+    // 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 = (Propagate->GetTec2l()*normal)(eZ);
+    double LGearProj = -(mTGear.Transposed() * vGroundNormal)(eZ);
 
     switch (eContactType) {
     case ctBOGEY:
-      // Project the height in the local coordinate frame of the strut to compute the actual compression
-      // length. The strut is assumed to be parallel to Z in the body frame.
-      compressLength = vGroundNormal(eZ) < 0.0 ? height / vGroundNormal(eZ) : 0.0;
+      compressLength = LGearProj > 0.0 ? height * normalZ / LGearProj : 0.0;
       break;
     case ctSTRUCTURE:
-      compressLength = -height;
+      compressLength = height * normalZ / DotProduct(normal, normal);
       break;
     }
 
@@ -300,71 +312,60 @@ FGColumnVector3& FGLGear::Force(void)
 
       WOW = true;
 
-      // [The next equation should really use the vector to the contact patch of
-      // the tire including the strut compression and not the original vWhlBodyVec.]
+      // The following equations use the vector to the tire contact patch
+      // including the strut compression.
+      FGColumnVector3 vWhlDisplVec;
+
+      switch(eContactType) {
+      case ctBOGEY:
+        vWhlDisplVec = mTGear * FGColumnVector3(0., 0., -compressLength);
+        break;
+      case ctSTRUCTURE:
+        vWhlDisplVec = compressLength * vGroundNormal;
+        break;
+      }
+
+      FGColumnVector3 vWhlContactVec = vWhlBodyVec + vWhlDisplVec;
+      vActingXYZn = vXYZn + Tb2s * vWhlDisplVec;
+      FGColumnVector3 vBodyWhlVel = Propagate->GetPQR() * vWhlContactVec;
+      vBodyWhlVel += Propagate->GetUVW() - Propagate->GetTec2b() * cvel;
 
-      FGColumnVector3 vWhlContactVec = vWhlBodyVec - FGColumnVector3(0., 0., compressLength);
-      vWhlVelVec      =  Propagate->GetPQR() * vWhlContactVec;
-      vWhlVelVec     +=  Propagate->GetUVW() - Propagate->GetTec2b() * cvel;
+      vWhlVelVec = mTGear.Transposed() * vBodyWhlVel;
 
       InitializeReporting();
       ComputeSteeringAngle();
       ComputeGroundCoordSys();
 
-      vLocalWhlVel = Tb2g * vWhlVelVec;
+      vLocalWhlVel = Transform().Transposed() * vBodyWhlVel;
 
-      compressSpeed = -vLocalWhlVel(eZ);
+      compressSpeed = -vLocalWhlVel(eX);
       if (eContactType == ctBOGEY)
-        // Project the compression speed in the local coordinate frame of the strut
-        compressSpeed /= -vGroundNormal(eZ);
+        compressSpeed /= LGearProj;
 
       ComputeVerticalStrutForce();
 
-      // Compute the forces in the wheel ground plane.
+      // Compute the friction coefficients in the wheel ground plane.
       if (eContactType == ctBOGEY) {
         ComputeSlipAngle();
         ComputeBrakeForceCoefficient();
         ComputeSideForceCoefficient();
-        double sign = vLocalWhlVel(eX)>0?1.0:(vLocalWhlVel(eX)<0?-1.0:0.0);
-        vLocalForce(eX) = - ((1.0 - TirePressureNorm) * 30 + vLocalForce(eZ) * BrakeFCoeff) * sign;
-        vLocalForce(eY) = vLocalForce(eZ) * FCoeff;
-      }
-      else if (eContactType == ctSTRUCTURE) {
-        FGColumnVector3 vSlipVec = vLocalWhlVel;
-        vSlipVec(eZ) = 0.;
-        vSlipVec.Normalize();
-        vLocalForce -= staticFCoeff * vLocalForce(eZ) * vSlipVec;
       }
 
-      // Lag and attenuate the XY-plane forces dependent on velocity. This code
-      // uses a lag filter, C/(s + C) where "C" is the filter coefficient. When
-      // "C" is chosen at the frame rate (in Hz), the jittering is significantly
-      // reduced. This is because the jitter is present *at* the execution rate.
-      // If a coefficient is set to something equal to or less than zero, the
-      // filter is bypassed.
-
-      if (LongForceLagFilterCoeff > 0) vLocalForce(eX) = LongForceFilter.execute(vLocalForce(eX));
-      if (LatForceLagFilterCoeff > 0)  vLocalForce(eY) = LatForceFilter.execute(vLocalForce(eY));
-
-      if ((fabs(vLocalWhlVel(eX)) <= RFRV) && RFRV > 0) vLocalForce(eX) *= fabs(vLocalWhlVel(eX))/RFRV;
-      if ((fabs(vLocalWhlVel(eY)) <= SFRV) && SFRV > 0) vLocalForce(eY) *= fabs(vLocalWhlVel(eY))/SFRV;
-
-      // End section for attenuating gear jitter
-
-      // Transform the forces back to the body frame and compute the moment.
-
-      vForce  = Tg2b * vLocalForce;
-      vMoment = vWhlContactVec * vForce;
+      // Prepare the Jacobians and the Lagrange multipliers for later friction
+      // forces calculations.
+      ComputeJacobian(vWhlContactVec);
 
     } else { // Gear is NOT compressed
 
       WOW = false;
       compressLength = 0.0;
       compressSpeed = 0.0;
+      WheelSlip = 0.0;
+      StrutForce = 0.0;
 
       // Let wheel spin down slowly
-      vLocalWhlVel(eX) -= 13.0*dT;
-      if (vLocalWhlVel(eX) < 0.0) vLocalWhlVel(eX) = 0.0;
+      vWhlVelVec(eX) -= 13.0*dT;
+      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;
@@ -373,61 +374,56 @@ FGColumnVector3& FGLGear::Force(void)
     }
   }
 
-  ReportTakeoffOrLanding();
+  if (!fdmex->GetTrimStatus()) {
+    ReportTakeoffOrLanding();
 
-  // Require both WOW and LastWOW to be true before checking crash conditions
-  // to allow the WOW flag to be used in terminating a scripted run.
-  if (WOW && lastWOW) CrashDetect();
+    // Require both WOW and LastWOW to be true before checking crash conditions
+    // to allow the WOW flag to be used in terminating a scripted run.
+    if (WOW && lastWOW) CrashDetect();
 
-  lastWOW = WOW;
+    lastWOW = WOW;
+  }
 
-  return vForce;
+  return FGForce::GetBodyForces();
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // Build a local "ground" coordinate system defined by
-//  eX : projection of the rolling direction on the ground
-//  eY : projection of the sliping direction on the ground
-//  eZ : normal to the ground
+//  eX : normal to the ground
+//  eY : projection of the rolling direction on the ground
+//  eZ : projection of the sliping direction on the ground
 
 void FGLGear::ComputeGroundCoordSys(void)
 {
-  if( fabs(vGroundNormal(eZ)) < 1e-3 ) {
-    cout << "BOMB: ";
-    cout << vGroundNormal(eX) << "/" << vGroundNormal(eY) << "/" << 
-         vGroundNormal(eZ) << " - " << SteerAngle << endl;
-    return;
-  }
-
-  // Compute the rolling direction projected on the ground
-  // It consists in finding a vector 'r' such that 'r' lies in the plane (w,z) and r.n = 0 (scalar
-  // product) where:
-  // 'n' is the normal to the ground,
-  // (x,y,z) are the directions defined in the body coord system
-  // and 'w' is 'x' rotated by the steering angle (SteerAngle) in the plane (x,y). 
-  // r = u * w + v * z and r.n = 0 => v/u = -w.n/z.n = a
-  // We also want u**2+v**2=1 and u > 0 (i.e. r orientated in the same 'direction' than w)
-  // after some arithmetic, one finds that :
-  double a = -(vGroundNormal(eX)*cos(SteerAngle)+vGroundNormal(eY)*sin(SteerAngle)) / vGroundNormal(eZ);
-  double u = 1. / sqrt(1. + a*a);
-  double v = a * u;
-  FGColumnVector3 vRollingGroundVec = FGColumnVector3(u * cos(SteerAngle), u * sin(SteerAngle), v);
-
-  // The sliping direction is the cross product multiplication of the ground normal and rolling
-  // directions
-  FGColumnVector3 vSlipGroundVec = vGroundNormal * vRollingGroundVec;
-
-  Tg2b(1,1) = vRollingGroundVec(eX);
-  Tg2b(2,1) = vRollingGroundVec(eY);
-  Tg2b(3,1) = vRollingGroundVec(eZ);
-  Tg2b(1,2) = vSlipGroundVec(eX);
-  Tg2b(2,2) = vSlipGroundVec(eY);
-  Tg2b(3,2) = vSlipGroundVec(eZ);
-  Tg2b(1,3) = vGroundNormal(eX);
-  Tg2b(2,3) = vGroundNormal(eY);
-  Tg2b(3,3) = vGroundNormal(eZ);
-
-  Tb2g = Tg2b.Transposed();
+  // Euler angles are built up to create a local frame to describe the forces
+  // applied to the gear by the ground. Here pitch, yaw and roll do not have
+  // 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)
+    vOrient(eYaw) = 0.;
+  else
+    vOrient(eYaw) = atan2(vGroundNormal(eY), vGroundNormal(eX));
+  
+  vOrient(eRoll) = 0.;
+  UpdateCustomTransformMatrix();
+
+  if (eContactType == ctBOGEY) {
+    // In the case of a bogey, the third angle "roll" is used to align the axis eY and eZ
+    // to the rolling and sliping direction respectively.
+    FGColumnVector3 updatedRollingAxis = Transform().Transposed() * mTGear
+                                       * FGColumnVector3(-sin(SteerAngle), cos(SteerAngle), 0.);
+
+    vOrient(eRoll) = atan2(updatedRollingAxis(eY), -updatedRollingAxis(eZ));
+    UpdateCustomTransformMatrix();
+  }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -439,7 +435,7 @@ void FGLGear::ComputeRetractionState(void)
     GearUp   = true;
     WOW      = false;
     GearDown = false;
-    vLocalWhlVel.InitMatrix();
+    vWhlVelVec.InitMatrix();
   } else if (gearPos > 0.99) {
     GearDown = true;
     GearUp   = false;
@@ -450,14 +446,13 @@ void FGLGear::ComputeRetractionState(void)
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+// Calculate tire slip angle.
 
 void FGLGear::ComputeSlipAngle(void)
 {
-  // Calculate tire slip angle.
-  WheelSlip = -atan2(vLocalWhlVel(eY), fabs(vLocalWhlVel(eX)))*radtodeg;
-
-  // Filter the wheel slip angle
-  if (WheelSlipLagFilterCoeff > 0) WheelSlip = WheelSlipFilter.execute(WheelSlip);
+// Check that the speed is non-null otherwise use the current angle
+  if (vLocalWhlVel.Magnitude(eY,eZ) > 1E-3)
+    WheelSlip = -atan2(vLocalWhlVel(eZ), fabs(vLocalWhlVel(eY)))*radtodeg;
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -474,7 +469,13 @@ void FGLGear::ComputeSteeringAngle(void)
     SteerAngle = 0.0;
     break;
   case stCaster:
-    SteerAngle = atan2(fabs(vWhlVelVec(eX)), vWhlVelVec(eY));
+    if (!Castered)
+      SteerAngle = degtorad * FCS->GetSteerPosDeg(GearNumber);
+    else {
+      // Check that the speed is non-null otherwise use the current angle
+      if (vWhlVelVec.Magnitude(eX,eY) > 0.1)
+        SteerAngle = atan2(vWhlVelVec(eY), fabs(vWhlVelVec(eX)));
+    }
     break;
   default:
     cerr << "Improper steering type membership detected for this gear." << endl;
@@ -529,20 +530,18 @@ void FGLGear::InitializeReporting(void)
 
 void FGLGear::ReportTakeoffOrLanding(void)
 {
-  double deltaT = State->Getdt()*Exec->GetGroundReactions()->GetRate();
-
   if (FirstContact)
-    LandingDistanceTraveled += Auxiliary->GetVground()*deltaT;
+    LandingDistanceTraveled += Auxiliary->GetVground()*dT;
 
   if (StartedGroundRun) {
-    TakeoffDistanceTraveled50ft += Auxiliary->GetVground()*deltaT;
-    if (WOW) TakeoffDistanceTraveled += Auxiliary->GetVground()*deltaT;
+    TakeoffDistanceTraveled50ft += Auxiliary->GetVground()*dT;
+    if (WOW) TakeoffDistanceTraveled += Auxiliary->GetVground()*dT;
   }
 
   if ( ReportEnable
        && Auxiliary->GetVground() <= 0.05
        && !LandingReported
-       && Exec->GetGroundReactions()->GetWOW())
+       && GroundReactions->GetWOW())
   {
     if (debug_lvl > 0) Report(erLand);
   }
@@ -550,7 +549,7 @@ void FGLGear::ReportTakeoffOrLanding(void)
   if ( ReportEnable
        && !TakeoffReported
        && (Propagate->GetDistanceAGL() - vLocalGear(eZ)) > 50.0
-       && !Exec->GetGroundReactions()->GetWOW())
+       && !GroundReactions->GetWOW())
   {
     if (debug_lvl > 0) Report(erTakeoff);
   }
@@ -564,12 +563,12 @@ void FGLGear::ReportTakeoffOrLanding(void)
 void FGLGear::CrashDetect(void)
 {
   if ( (compressLength > 500.0 ||
-      vForce.Magnitude() > 100000000.0 ||
-      vMoment.Magnitude() > 5000000000.0 ||
-      SinkRate > 1.4666*30 ) && !State->IntegrationSuspended())
+      vFn.Magnitude() > 100000000.0 ||
+      GetMoments().Magnitude() > 5000000000.0 ||
+      SinkRate > 1.4666*30 ) && !fdmex->IntegrationSuspended())
   {
     PutMessage("Crash Detected: Simulation FREEZE.");
-    State->SuspendIntegration();
+    fdmex->SuspendIntegration();
   }
 }
 
@@ -667,10 +666,10 @@ void FGLGear::ComputeVerticalStrutForce(void)
   switch (eContactType) {
   case ctBOGEY:
     // Project back the strut force in the local coordinate frame of the ground
-    vLocalForce(eZ) =  StrutForce / vGroundNormal(eZ);
+    vFn(eX) = StrutForce / (mTGear.Transposed()*vGroundNormal)(eZ);
     break;
   case ctSTRUCTURE:
-    vLocalForce(eZ) = -StrutForce;
+    vFn(eX) = -StrutForce;
     break;
   }
 
@@ -691,6 +690,99 @@ double FGLGear::GetGearUnitPos(void)
   return GearPos;
 }
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+// Compute the jacobian entries for the friction forces resolution later
+// in FGPropagate
+
+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)) {
+    FGColumnVector3 velocityDirection = vLocalWhlVel;
+
+    StaticFriction = false;
+
+    velocityDirection(eX) = 0.;
+    velocityDirection.Normalize();
+
+    LMultiplier[ftDynamic].ForceJacobian = Transform()*velocityDirection;
+    LMultiplier[ftDynamic].MomentJacobian = vWhlContactVec * LMultiplier[ftDynamic].ForceJacobian;
+    LMultiplier[ftDynamic].Max = 0.;
+    LMultiplier[ftDynamic].Min = -fabs(dynamicFCoeff * vFn(eX));
+    LMultiplier[ftDynamic].value = Constrain(LMultiplier[ftDynamic].Min, LMultiplier[ftDynamic].value, LMultiplier[ftDynamic].Max);
+  }
+  else {
+    // Static friction is used for ctSTRUCTURE when the contact point is not moving.
+    // It is always used for ctBOGEY elements because the friction coefficients
+    // of a tyre depend on the direction of the movement (roll & side directions).
+    // This cannot be handled properly by the so-called "dynamic friction".
+    StaticFriction = true;
+
+    LMultiplier[ftRoll].ForceJacobian = Transform()*FGColumnVector3(0.,1.,0.);
+    LMultiplier[ftSide].ForceJacobian = Transform()*FGColumnVector3(0.,0.,1.);
+    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[ftSide].Max = fabs(FCoeff * vFn(eX));
+      break;
+    case ctSTRUCTURE:
+      LMultiplier[ftRoll].Max = fabs(staticFCoeff * vFn(eX));
+      LMultiplier[ftSide].Max = fabs(staticFCoeff * vFn(eX));
+      break;
+    }
+
+    LMultiplier[ftRoll].Min = -LMultiplier[ftRoll].Max;
+    LMultiplier[ftSide].Min = -LMultiplier[ftSide].Max;
+    LMultiplier[ftRoll].value = Constrain(LMultiplier[ftRoll].Min, LMultiplier[ftRoll].value, LMultiplier[ftRoll].Max);
+    LMultiplier[ftSide].value = Constrain(LMultiplier[ftSide].Min, LMultiplier[ftSide].value, LMultiplier[ftSide].Max);
+  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+// This function is used by the MultiplierIterator class to enumerate the
+// Lagrange multipliers of a landing gear. This allows to encapsulate the storage
+// of the multipliers in FGLGear without exposing it. From an outside point of
+// view, each FGLGear instance has a number of Lagrange multipliers which can be
+// accessed through this routine without knowing the exact constraint which they
+// model.
+
+FGPropagate::LagrangeMultiplier* FGLGear::GetMultiplierEntry(int entry)
+{
+  switch(entry) {
+  case 0:
+    if (StaticFriction)
+      return &LMultiplier[ftRoll];
+    else
+      return &LMultiplier[ftDynamic];
+  case 1:
+    if (StaticFriction)
+      return &LMultiplier[ftSide];
+  default:
+    return NULL;
+  }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+// This routine is called after the Lagrange multiplier has been computed. The
+// friction forces of the landing gear are then updated accordingly.
+FGColumnVector3& FGLGear::UpdateForces(void)
+{
+  if (StaticFriction) {
+    vFn(eY) = LMultiplier[ftRoll].value;
+    vFn(eZ) = LMultiplier[ftSide].value;
+  }
+  else
+    vFn += LMultiplier[ftDynamic].value * (Transform ().Transposed() * LMultiplier[ftDynamic].ForceJacobian);
+
+  // Return the updated force in the body frame
+  return FGForce::GetBodyForces();
+}
+
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 void FGLGear::bind(void)
@@ -700,30 +792,39 @@ void FGLGear::bind(void)
   base_property_name = CreateIndexedPropertyName("gear/unit", GearNumber);
   if (eContactType == ctBOGEY) {
     property_name = base_property_name + "/slip-angle-deg";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), &WheelSlip );
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &WheelSlip );
     property_name = base_property_name + "/WOW";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), &WOW );
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &WOW );
     property_name = base_property_name + "/wheel-speed-fps";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), (FGLGear*)this,
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), (FGLGear*)this,
                           &FGLGear::GetWheelRollVel);
     property_name = base_property_name + "/z-position";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), (FGLGear*)this,
-                          &FGLGear::GetZPosition, &FGLGear::SetZPosition);
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), (FGForce*)this,
+                          &FGForce::GetLocationZ, &FGForce::SetLocationZ);
     property_name = base_property_name + "/compression-ft";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), &compressLength );
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &compressLength );
     property_name = base_property_name + "/side_friction_coeff";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), &FCoeff );
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &FCoeff );
 
     property_name = base_property_name + "/static_friction_coeff";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), &staticFCoeff );
-
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &staticFCoeff );
+    property_name = base_property_name + "/rolling_friction_coeff";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &rollingFCoeff );
+    property_name = base_property_name + "/dynamic_friction_coeff";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &dynamicFCoeff );
+
+    if (eSteerType == stCaster) {
+      property_name = base_property_name + "/steering-angle-deg";
+      fdmex->GetPropertyManager()->Tie( property_name.c_str(), this, &FGLGear::GetSteerAngleDeg );
+      property_name = base_property_name + "/castered";
+      fdmex->GetPropertyManager()->Tie( property_name.c_str(), &Castered);
+    }
   }
 
   if( isRetractable ) {
     property_name = base_property_name + "/pos-norm";
-    Exec->GetPropertyManager()->Tie( property_name.c_str(), &GearPos );
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &GearPos );
   }
-
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -735,7 +836,7 @@ void FGLGear::Report(ReportType repType)
   switch(repType) {
   case erLand:
     cout << endl << "Touchdown report for " << name << " (WOW at time: "
-         << Exec->GetState()->Getsim_time() << " seconds)" << endl;
+         << fdmex->GetSimTime() << " seconds)" << endl;
     cout << "  Sink rate at contact:  " << SinkRate                << " fps,    "
                                 << SinkRate*0.3048          << " mps"     << endl;
     cout << "  Contact ground speed:  " << GroundSpeed*.5925       << " knots,  "
@@ -750,18 +851,20 @@ void FGLGear::Report(ReportType repType)
     break;
   case erTakeoff:
     cout << endl << "Takeoff report for " << name << " (Liftoff at time: "
-         << Exec->GetState()->Getsim_time() << " seconds)" << endl;
+        << fdmex->GetSimTime() << " seconds)" << endl;
     cout << "  Distance traveled:                " << TakeoffDistanceTraveled
          << " ft,     " << TakeoffDistanceTraveled*0.3048  << " meters"  << endl;
     cout << "  Distance traveled (over 50'):     " << TakeoffDistanceTraveled50ft
          << " ft,     " << TakeoffDistanceTraveled50ft*0.3048 << " meters" << endl;
-    cout << "  [Altitude (ASL): " << Exec->GetPropagate()->GetAltitudeASL() << " ft. / "
-         << Exec->GetPropagate()->GetAltitudeASLmeters() << " m  | Temperature: "
-         << Exec->GetAtmosphere()->GetTemperature() - 459.67 << " F / "
-         << RankineToCelsius(Exec->GetAtmosphere()->GetTemperature()) << " C]" << endl;
-    cout << "  [Velocity (KCAS): " << Exec->GetAuxiliary()->GetVcalibratedKTS() << "]" << endl;
+    cout << "  [Altitude (ASL): " << Propagate->GetAltitudeASL() << " ft. / "
+         << Propagate->GetAltitudeASLmeters() << " m  | Temperature: "
+         << fdmex->GetAtmosphere()->GetTemperature() - 459.67 << " F / "
+         << RankineToCelsius(fdmex->GetAtmosphere()->GetTemperature()) << " C]" << endl;
+    cout << "  [Velocity (KCAS): " << Auxiliary->GetVcalibratedKTS() << "]" << endl;
     TakeoffReported = true;
     break;
+  case erNone:
+    break;
   }
 }
 
@@ -791,7 +894,7 @@ void FGLGear::Debug(int from)
   if (debug_lvl & 1) { // Standard console startup message output
     if (from == 0) { // Constructor - loading and initialization
       cout << "    " << sContactType << " " << name          << endl;
-      cout << "      Location: "         << vXYZ          << endl;
+      cout << "      Location: "         << vXYZn          << endl;
       cout << "      Spring Constant:  " << kSpring       << endl;
 
       if (eDampType == dtLinear)
@@ -812,9 +915,6 @@ void FGLGear::Debug(int from)
         cout << "      Grouping:         " << sBrakeGroup   << endl;
         cout << "      Max Steer Angle:  " << maxSteerAngle << endl;
         cout << "      Retractable:      " << isRetractable  << endl;
-        cout << "      Relaxation Velocities:" << endl;
-        cout << "        Rolling:          " << RFRV << endl;
-        cout << "        Side:             " << SFRV << endl;
       }
     }
   }