]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/JSBSim/models/FGLGear.cpp
Merge branch 'jmt/units-fix' into maint
[flightgear.git] / src / FDM / JSBSim / models / FGLGear.cpp
index 08011197629defd73c05bcd6c18e418a22b64d68..b0c448b41777baf98a5ff7375f6017505f067c27 100644 (file)
@@ -3,11 +3,12 @@
  Module:       FGLGear.cpp
  Author:       Jon S. Berndt
                Norman H. Princen
+               Bertrand Coconnier
  Date started: 11/18/99
  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
@@ -33,12 +34,23 @@ HISTORY
 --------------------------------------------------------------------------------
 11/18/99   JSB   Created
 01/30/01   NHP   Extended gear model to properly simulate steering and braking
+07/08/09   BC    Modified gear model to support large angles between aircraft and ground
 
 /%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 INCLUDES
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 #include "FGLGear.h"
+#include "FGState.h"
+#include "FGGroundReactions.h"
+#include "FGFCS.h"
+#include "FGAuxiliary.h"
+#include "FGAtmosphere.h"
+#include "FGMassBalance.h"
+#include "math/FGTable.h"
+#include <cstdlib>
+
+using namespace std;
 
 namespace JSBSim {
 
@@ -53,13 +65,18 @@ GLOBAL DATA
 static const char *IdSrc = "$Id$";
 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)
 {
   Element *force_table=0;
   Element *dampCoeff=0;
@@ -79,7 +96,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"))
@@ -123,7 +141,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;
     }
@@ -137,8 +155,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;
@@ -176,7 +223,13 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
     }
   }
 
-  State = Exec->GetState();
+  State       = fdmex->GetState();
+  Aircraft    = fdmex->GetAircraft();
+  Propagate   = fdmex->GetPropagate();
+  Auxiliary   = fdmex->GetAuxiliary();
+  FCS         = fdmex->GetFCS();
+  MassBalance = fdmex->GetMassBalance();
+
   LongForceLagFilterCoeff = 1/State->Getdt(); // default longitudinal force filter coefficient
   LatForceLagFilterCoeff  = 1/State->Getdt(); // default lateral force filter coefficient
 
@@ -190,6 +243,9 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
     }
   }
 
+  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");
@@ -197,6 +253,8 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
     WheelSlipLagFilterCoeff = wheel_slip_angle_lag_elem->GetDataAsNumber();
   }
   
+  WheelSlipFilter = Filter(WheelSlipLagFilterCoeff, State->Getdt());
+
   GearUp = false;
   GearDown = true;
   GearPos  = 1.0;
@@ -206,13 +264,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;
@@ -220,12 +271,11 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
   TakeoffReported = LandingReported = false;
   LandingDistanceTraveled = TakeoffDistanceTraveled = TakeoffDistanceTraveled50ft = 0.0;
   MaximumStrutForce = MaximumStrutTravel = 0.0;
-  SideForce = RollingForce = 0.0;
   SinkRate = GroundSpeed = 0.0;
 
-  vWhlBodyVec = MassBalance->StructuralToBody(vXYZ);
-
+  vWhlBodyVec = MassBalance->StructuralToBody(vXYZn);
   vLocalGear = Propagate->GetTb2l() * vWhlBodyVec;
+  vWhlVelVec.InitMatrix();
 
   compressLength  = 0.0;
   compressSpeed   = 0.0;
@@ -235,122 +285,60 @@ FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number) :
   WheelSlip = 0.0;
   TirePressureNorm = 1.0;
 
-  SideWhlVel    = 0.0;
-  RollingWhlVel = 0.0;
-
-  SinWheel = 0.0;
-  CosWheel = 0.0;
+  // Set Pacejka terms
 
-  prevSlipIn  = 0.0;
-  prevSlipOut = 0.0;
+  Stiffness = 0.06;
+  Shape = 2.8;
+  Peak = staticFCoeff;
+  Curvature = 1.03;
 
   Debug(0);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-FGLGear::FGLGear(const FGLGear& lgear)
-{
-  GearNumber = lgear.GearNumber;
-  State    = lgear.State;
-  Aircraft = lgear.Aircraft;
-  Propagate = lgear.Propagate;
-  Auxiliary = lgear.Auxiliary;
-  Exec     = lgear.Exec;
-  FCS      = lgear.FCS;
-  MassBalance = lgear.MassBalance;
-
-  vXYZ = lgear.vXYZ;
-  vMoment = lgear.vMoment;
-  vWhlBodyVec = lgear.vWhlBodyVec;
-  vLocalGear = lgear.vLocalGear;
-
-  WOW                = lgear.WOW;
-  lastWOW            = lgear.lastWOW;
-  ReportEnable       = lgear.ReportEnable;
-  FirstContact       = lgear.FirstContact;
-  StartedGroundRun   = lgear.StartedGroundRun;
-  LandingDistanceTraveled   = lgear.LandingDistanceTraveled;
-  TakeoffDistanceTraveled   = lgear.TakeoffDistanceTraveled;
-  TakeoffDistanceTraveled50ft   = lgear.TakeoffDistanceTraveled50ft;
-  MaximumStrutForce  = lgear.MaximumStrutForce;
-  MaximumStrutTravel = lgear.MaximumStrutTravel;
-  SideForce          = lgear.SideForce;
-  RollingForce       = lgear.RollingForce;
-
-  kSpring         = lgear.kSpring;
-  bDamp           = lgear.bDamp;
-  bDampRebound    = lgear.bDampRebound;
-  compressLength  = lgear.compressLength;
-  compressSpeed   = lgear.compressSpeed;
-  staticFCoeff    = lgear.staticFCoeff;
-  dynamicFCoeff   = lgear.dynamicFCoeff;
-  rollingFCoeff   = lgear.rollingFCoeff;
-  brakePct        = lgear.brakePct;
-  maxCompLen      = lgear.maxCompLen;
-  SinkRate        = lgear.SinkRate;
-  GroundSpeed     = lgear.GroundSpeed;
-  LandingReported = lgear.LandingReported;
-  TakeoffReported = lgear.TakeoffReported;
-  name            = lgear.name;
-  sSteerType      = lgear.sSteerType;
-  sRetractable    = lgear.sRetractable;
-  sContactType    = lgear.sContactType;
-  eContactType    = lgear.eContactType;
-  sBrakeGroup     = lgear.sBrakeGroup;
-  eSteerType      = lgear.eSteerType;
-  eBrakeGrp       = lgear.eBrakeGrp;
-  maxSteerAngle   = lgear.maxSteerAngle;
-  isRetractable   = lgear.isRetractable;
-  GearUp          = lgear.GearUp;
-  GearDown        = lgear.GearDown;
-  GearPos         = lgear.GearPos;
-  useFCSGearPos   = lgear.useFCSGearPos;
-  WheelSlip       = lgear.WheelSlip;
-  TirePressureNorm = lgear.TirePressureNorm;
-  Servicable      = lgear.Servicable;
-  ForceY_Table    = lgear.ForceY_Table;
-  CosWheel        = lgear.CosWheel;
-  SinWheel        = lgear.SinWheel;
-  prevOut         = lgear.prevOut;
-  prevIn          = lgear.prevIn;
-  prevSlipIn      = lgear.prevSlipIn;
-  prevSlipOut     = lgear.prevSlipOut;
-  RFRV            = lgear.RFRV;
-  SFRV            = lgear.SFRV;
-  LongForceLagFilterCoeff = lgear.LongForceLagFilterCoeff;
-  LatForceLagFilterCoeff = lgear.LatForceLagFilterCoeff;
-  WheelSlipLagFilterCoeff = lgear.WheelSlipLagFilterCoeff;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 FGLGear::~FGLGear()
 {
+  delete ForceY_Table;
   Debug(1);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-FGColumnVector3& FGLGear::Force(void)
+FGColumnVector3& FGLGear::GetBodyForces(void)
 {
-  double t = Exec->GetState()->Getsim_time();
-  dT = State->Getdt()*Exec->GetGroundReactions()->GetRate();
+  double t = fdmex->GetState()->Getsim_time();
+  dT = State->Getdt()*fdmex->GetGroundReactions()->GetRate();
 
-  vForce.InitMatrix();
-  vMoment.InitMatrix();
+  vFn.InitMatrix();
 
   if (isRetractable) ComputeRetractionState();
 
   if (GearDown) {
+    double verticalZProj = 0.;
 
-    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);
-    compressLength = -Exec->GetGroundCallback()->GetAGLevel(t, gearLoc, contact, normal, cvel);
-
-    // The compression length is measured in the Z-axis, only, at this time.
+    // 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);
+    vGroundNormal = -1. * Propagate->GetTec2b() * normal;
+
+    // The height returned above 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)
+    // of in the normal direction to the ground (STRUCTURE)
+    switch (eContactType) {
+    case ctBOGEY:
+      verticalZProj = (Propagate->GetTb2l()*mTGear*FGColumnVector3(0.,0.,1.))(eZ);
+      compressLength = verticalZProj > 0.0 ? -height / verticalZProj : 0.0;
+      break;
+    case ctSTRUCTURE:
+      verticalZProj = (Propagate->GetTec2l()*normal)(eZ);
+      compressLength = fabs(verticalZProj) > 0.0 ? -height / verticalZProj : 0.0;
+      break;
+    }
 
     if (compressLength > 0.00) {
 
@@ -359,73 +347,73 @@ FGColumnVector3& FGLGear::Force(void)
       // [The next equation should really use the vector to the contact patch of
       // the tire including the strut compression and not the original vWhlBodyVec.]
 
-      vWhlVelVec      =  Propagate->GetTb2l() * (Propagate->GetPQR() * vWhlBodyVec);
-      vWhlVelVec     +=  Propagate->GetVel() - cvel;
-      compressSpeed   =  vWhlVelVec(eZ);
+      FGColumnVector3 vWhlDisplVec = mTGear * FGColumnVector3(0., 0., compressLength);
+      FGColumnVector3 vWhlContactVec = vWhlBodyVec - vWhlDisplVec;
+      vActingXYZn = vXYZn - Tb2s * vWhlDisplVec;
+      FGColumnVector3 vBodyWhlVel  = Propagate->GetPQR() * vWhlContactVec;
+      vBodyWhlVel += Propagate->GetUVW() - Propagate->GetTec2b() * cvel;
+
+      vWhlVelVec = mTGear.Transposed() * vBodyWhlVel;
 
       InitializeReporting();
-      ComputeBrakeForceCoefficient();
       ComputeSteeringAngle();
-      ComputeSlipAngle();
-      ComputeSideForceCoefficient();
-      ComputeVerticalStrutForce();
+      ComputeGroundCoordSys();
 
-      // Compute the forces in the wheel ground plane.
-
-      double sign = RollingWhlVel>0?1.0:(RollingWhlVel<0?-1.0:0.0);
-      RollingForce = ((1.0 - TirePressureNorm) * 30 + vLocalForce(eZ) * BrakeFCoeff) * sign;
-      SideForce    = vLocalForce(eZ) * FCoeff;
-
-      // Transform these forces back to the local reference frame.
-
-      vLocalForce(eX) = RollingForce*CosWheel - SideForce*SinWheel;
-      vLocalForce(eY) = SideForce*CosWheel    + RollingForce*SinWheel;
-
-      // Transform the forces back to the body frame and compute the moment.
+      vLocalWhlVel = Transform().Transposed() * vBodyWhlVel;
 
-      vForce  = Propagate->GetTl2b() * vLocalForce;
-
-      // Lag and attenuate the XY-plane forces dependent on velocity
-
-      double ca, cb, denom;
-      FGColumnVector3 Output;
+      switch (eContactType) {
+      case ctBOGEY:
+        // Compression speed along the strut
+        compressSpeed = -vWhlVelVec(eZ);
+      case ctSTRUCTURE:
+        // Compression speed along the ground normal
+        compressSpeed = -vLocalWhlVel(eX);
+      }
 
-      // This code implements a lag filter, C/(s + C) where
-      // "C" is the filter coefficient. When "C" is chosen at the 
-      // frame rate (in Hz), the jittering is significantly reduced. This is because
-      // the jitter is present *at* the execution rate.
-      // If a coefficient is set to something equal to or less than zero, the filter
-      // is bypassed.
+      ComputeVerticalStrutForce();
 
-      if (LongForceLagFilterCoeff > 0) { 
-        denom = 2.00 + dT*LongForceLagFilterCoeff;
-        ca = dT*LongForceLagFilterCoeff / denom;
-        cb = (2.00 - dT*LongForceLagFilterCoeff) / denom;
-        Output(eX) = vForce(eX) * ca + prevIn(eX) * ca + prevOut(eX) * cb;
-        vForce(eX) = Output(eX);
+      // Compute the forces in the wheel ground plane.
+      if (eContactType == ctBOGEY) {
+        ComputeSlipAngle();
+        ComputeBrakeForceCoefficient();
+        ComputeSideForceCoefficient();
+        double sign = vLocalWhlVel(eY)>0?1.0:(vLocalWhlVel(eY)<0?-1.0:0.0);
+        vFn(eY) = - ((1.0 - TirePressureNorm) * 30 + vFn(eX) * BrakeFCoeff) * sign;
+        vFn(eZ) = vFn(eX) * FCoeff;
       }
-      if (LatForceLagFilterCoeff > 0) { 
-        denom = 2.00 + dT*LatForceLagFilterCoeff;
-        ca = dT*LatForceLagFilterCoeff / denom;
-        cb = (2.00 - dT*LatForceLagFilterCoeff) / denom;
-        Output(eY) = vForce(eY) * ca + prevIn(eY) * ca + prevOut(eY) * cb;
-        vForce(eY) = Output(eY);
+      else if (eContactType == ctSTRUCTURE) {
+        FGColumnVector3 vSlipVec = vLocalWhlVel;
+        vSlipVec(eX) = 0.;
+        vSlipVec.Normalize();
+        vFn -= staticFCoeff * vFn(eX) * vSlipVec;
       }
 
-      prevIn = vForce;
-      prevOut = Output;
+      // 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 ((fabs(RollingWhlVel) <= RFRV) && RFRV > 0) vForce(eX) *= fabs(RollingWhlVel)/RFRV;
-      if ((fabs(SideWhlVel) <= SFRV) && SFRV > 0) vForce(eY) *= fabs(SideWhlVel)/SFRV;
+      if (LongForceLagFilterCoeff > 0) vFn(eY) = LongForceFilter.execute(vFn(eY));
+      if (LatForceLagFilterCoeff > 0)  vFn(eZ) = LatForceFilter.execute(vFn(eZ));
 
-  // End section for attentuating gear jitter
+      if ((fabs(vLocalWhlVel(eY)) <= RFRV) && RFRV > 0) vFn(eY) *= fabs(vLocalWhlVel(eY))/RFRV;
+      if ((fabs(vLocalWhlVel(eZ)) <= SFRV) && SFRV > 0) vFn(eZ) *= fabs(vLocalWhlVel(eZ))/SFRV;
 
-      vMoment = vWhlBodyVec * vForce;
+      // End section for attenuating gear jitter
 
     } else { // Gear is NOT compressed
 
       WOW = false;
       compressLength = 0.0;
+      compressSpeed = 0.0;
+      WheelSlip = 0.0;
+      StrutForce = 0.0;
+
+      // Let wheel spin down slowly
+      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;
@@ -442,7 +430,46 @@ FGColumnVector3& FGLGear::Force(void)
 
   lastWOW = WOW;
 
-  return vForce;
+  return FGForce::GetBodyForces();
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+// Build a local "ground" coordinate system defined by
+//  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)
+{
+  // 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();
+  }
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -454,6 +481,7 @@ void FGLGear::ComputeRetractionState(void)
     GearUp   = true;
     WOW      = false;
     GearDown = false;
+    vWhlVelVec.InitMatrix();
   } else if (gearPos > 0.99) {
     GearDown = true;
     GearUp   = false;
@@ -467,27 +495,11 @@ void FGLGear::ComputeRetractionState(void)
 
 void FGLGear::ComputeSlipAngle(void)
 {
-  // Transform the wheel velocities from the local axis system to the wheel axis system.
-  RollingWhlVel = vWhlVelVec(eX)*CosWheel + vWhlVelVec(eY)*SinWheel;
-  SideWhlVel    = vWhlVelVec(eY)*CosWheel - vWhlVelVec(eX)*SinWheel;
-
   // Calculate tire slip angle.
-    WheelSlip = atan2(SideWhlVel, fabs(RollingWhlVel))*radtodeg;
+  WheelSlip = -atan2(vLocalWhlVel(eZ), fabs(vLocalWhlVel(eY)))*radtodeg;
 
   // Filter the wheel slip angle
-
-  double SlipOutput, ca, cb, denom;
-
-  if (WheelSlipLagFilterCoeff > 0) {
-    denom = 2.00 + dT*WheelSlipLagFilterCoeff;
-    ca = dT*WheelSlipLagFilterCoeff / denom;
-    cb = (2.00 - dT*WheelSlipLagFilterCoeff) / denom;
-
-    SlipOutput = ca * (WheelSlip + prevSlipIn) + cb * prevSlipOut;
-
-    prevSlipIn = WheelSlip;
-    WheelSlip = prevSlipOut = SlipOutput;
-  }
+  if (WheelSlipLagFilterCoeff > 0) WheelSlip = WheelSlipFilter.execute(WheelSlip);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -496,9 +508,6 @@ void FGLGear::ComputeSlipAngle(void)
 
 void FGLGear::ComputeSteeringAngle(void)
 {
-  double casterLocalFrameAngleRad = 0.0;
-  double casterAngle = 0.0;
-
   switch (eSteerType) {
   case stSteer:
     SteerAngle = degtorad * FCS->GetSteerPosDeg(GearNumber);
@@ -507,20 +516,12 @@ void FGLGear::ComputeSteeringAngle(void)
     SteerAngle = 0.0;
     break;
   case stCaster:
-    // This is not correct for castering gear. Should make steer angle parallel
-    // to the actual velocity vector of the wheel, given aircraft velocity vector
-    // and omega.
-    SteerAngle = 0.0;
-    casterLocalFrameAngleRad = acos(vWhlVelVec(eX)/vWhlVelVec.Magnitude());
-    casterAngle = casterLocalFrameAngleRad - Propagate->GetEuler(ePsi);
+    SteerAngle = atan2(vWhlVelVec(eY), fabs(vWhlVelVec(eX)));
     break;
   default:
     cerr << "Improper steering type membership detected for this gear." << endl;
     break;
   }
-
-  SinWheel      = sin(Propagate->GetEuler(ePsi) + SteerAngle);
-  CosWheel      = cos(Propagate->GetEuler(ePsi) + SteerAngle);
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -570,7 +571,7 @@ void FGLGear::InitializeReporting(void)
 
 void FGLGear::ReportTakeoffOrLanding(void)
 {
-  double deltaT = State->Getdt()*Exec->GetGroundReactions()->GetRate();
+  double deltaT = State->Getdt()*fdmex->GetGroundReactions()->GetRate();
 
   if (FirstContact)
     LandingDistanceTraveled += Auxiliary->GetVground()*deltaT;
@@ -583,7 +584,7 @@ void FGLGear::ReportTakeoffOrLanding(void)
   if ( ReportEnable
        && Auxiliary->GetVground() <= 0.05
        && !LandingReported
-       && Exec->GetGroundReactions()->GetWOW())
+       && fdmex->GetGroundReactions()->GetWOW())
   {
     if (debug_lvl > 0) Report(erLand);
   }
@@ -591,7 +592,7 @@ void FGLGear::ReportTakeoffOrLanding(void)
   if ( ReportEnable
        && !TakeoffReported
        && (Propagate->GetDistanceAGL() - vLocalGear(eZ)) > 50.0
-       && !Exec->GetGroundReactions()->GetWOW())
+       && !fdmex->GetGroundReactions()->GetWOW())
   {
     if (debug_lvl > 0) Report(erTakeoff);
   }
@@ -605,8 +606,8 @@ void FGLGear::ReportTakeoffOrLanding(void)
 void FGLGear::CrashDetect(void)
 {
   if ( (compressLength > 500.0 ||
-      vForce.Magnitude() > 100000000.0 ||
-      vMoment.Magnitude() > 5000000000.0 ||
+      vFn.Magnitude() > 100000000.0 ||
+      GetMoments().Magnitude() > 5000000000.0 ||
       SinkRate > 1.4666*30 ) && !State->IntegrationSuspended())
   {
     PutMessage("Crash Detected: Simulation FREEZE.");
@@ -655,27 +656,22 @@ void FGLGear::ComputeBrakeForceCoefficient(void)
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-// Compute the sideforce coefficients using similar assumptions to LaRCSim for now.
-// Allow a maximum of 10 degrees tire slip angle before wheel slides.  At that point,
-// transition from static to dynamic friction.  There are more complicated formulations
-// of this that avoid the discrete jump (similar to Pacejka).  Will fix this later.
+// Compute the sideforce coefficients using Pacejka's Magic Formula.
+//
+//   y(x) = D sin {C arctan [Bx - E(Bx - arctan Bx)]}
+//
+// Where: B = Stiffness Factor (0.06, here)
+//        C = Shape Factor (2.8, here)
+//        D = Peak Factor (0.8, here)
+//        E = Curvature Factor (1.03, here)
 
 void FGLGear::ComputeSideForceCoefficient(void)
 {
   if (ForceY_Table) {
-
     FCoeff = ForceY_Table->GetValue(WheelSlip);
-
   } else {
-
-    if (fabs(WheelSlip) <= 10.0) {
-      FCoeff = staticFCoeff*WheelSlip/10.0;
-    } else if (fabs(WheelSlip) <= 40.0) {
-      FCoeff = (dynamicFCoeff*(fabs(WheelSlip) - 10.0)/10.0
-                + staticFCoeff*(40.0 - fabs(WheelSlip))/10.0)*(WheelSlip>=0?1.0:-1.0);
-    } else {
-      FCoeff = dynamicFCoeff*(WheelSlip>=0?1.0:-1.0);
-    }
+    double StiffSlip = Stiffness*WheelSlip;
+    FCoeff = Peak * sin(Shape*atan(StiffSlip - Curvature*(StiffSlip - atan(StiffSlip))));
   }
 }
 
@@ -706,10 +702,22 @@ void FGLGear::ComputeVerticalStrutForce(void)
       dampForce   =  compressSpeed * compressSpeed * bDampRebound;
 
   }
-  vLocalForce(eZ) =  min(springForce + dampForce, (double)0.0);
+
+  StrutForce = min(springForce + dampForce, (double)0.0);
+
+  // The reaction force of the wheel is always normal to the ground
+  switch (eContactType) {
+  case ctBOGEY:
+    // Project back the strut force in the local coordinate frame of the ground
+    vFn(eX) = StrutForce / (mTGear.Transposed()*vGroundNormal)(eZ);
+    break;
+  case ctSTRUCTURE:
+    vFn(eX) = -StrutForce;
+    break;
+  }
 
   // Remember these values for reporting
-  MaximumStrutForce = max(MaximumStrutForce, fabs(vLocalForce(eZ)));
+  MaximumStrutForce = max(MaximumStrutForce, fabs(StrutForce));
   MaximumStrutTravel = max(MaximumStrutTravel, fabs(compressLength));
 }
 
@@ -729,26 +737,38 @@ double FGLGear::GetGearUnitPos(void)
 
 void FGLGear::bind(void)
 {
-  char property_name[80];
+  string property_name;
+  string base_property_name;
+  base_property_name = CreateIndexedPropertyName("gear/unit", GearNumber);
   if (eContactType == ctBOGEY) {
-    snprintf(property_name, 80, "gear/unit[%d]/slip-angle-deg", GearNumber);
-    Exec->GetPropertyManager()->Tie( property_name, &WheelSlip );
-    snprintf(property_name, 80, "gear/unit[%d]/WOW", GearNumber);
-    Exec->GetPropertyManager()->Tie( property_name, &WOW );
-    snprintf(property_name, 80, "gear/unit[%d]/wheel-speed-fps", GearNumber);
-    Exec->GetPropertyManager()->Tie( property_name, &RollingWhlVel );
-    snprintf(property_name, 80, "gear/unit[%d]/z-position", GearNumber);
-    Exec->GetPropertyManager()->Tie( property_name, (FGLGear*)this,
-                          &FGLGear::GetZPosition, &FGLGear::SetZPosition);
-    snprintf(property_name, 80, "gear/unit[%d]/compression-ft", GearNumber);
-    Exec->GetPropertyManager()->Tie( property_name, &compressLength );
+    property_name = base_property_name + "/slip-angle-deg";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &WheelSlip );
+    property_name = base_property_name + "/WOW";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &WOW );
+    property_name = base_property_name + "/wheel-speed-fps";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), (FGLGear*)this,
+                          &FGLGear::GetWheelRollVel);
+    property_name = base_property_name + "/z-position";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), (FGForce*)this,
+                          &FGForce::GetLocationZ, &FGForce::SetLocationZ);
+    property_name = base_property_name + "/compression-ft";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &compressLength );
+    property_name = base_property_name + "/side_friction_coeff";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &FCoeff );
+
+    property_name = base_property_name + "/static_friction_coeff";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &staticFCoeff );
+
+    if (eSteerType == stCaster) {
+      property_name = base_property_name + "/steering-angle-rad";
+      fdmex->GetPropertyManager()->Tie( property_name.c_str(), &SteerAngle );
+    }
   }
 
   if( isRetractable ) {
-    snprintf(property_name, 80, "gear/unit[%d]/pos-norm", GearNumber);
-    Exec->GetPropertyManager()->Tie( property_name, &GearPos );
+    property_name = base_property_name + "/pos-norm";
+    fdmex->GetPropertyManager()->Tie( property_name.c_str(), &GearPos );
   }
-
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -759,7 +779,8 @@ void FGLGear::Report(ReportType repType)
 
   switch(repType) {
   case erLand:
-    cout << endl << "Touchdown report for " << name << endl;
+    cout << endl << "Touchdown report for " << name << " (WOW at time: "
+         << fdmex->GetState()->Getsim_time() << " seconds)" << endl;
     cout << "  Sink rate at contact:  " << SinkRate                << " fps,    "
                                 << SinkRate*0.3048          << " mps"     << endl;
     cout << "  Contact ground speed:  " << GroundSpeed*.5925       << " knots,  "
@@ -773,13 +794,21 @@ void FGLGear::Report(ReportType repType)
     LandingReported = true;
     break;
   case erTakeoff:
-    cout << endl << "Takeoff report for " << name << endl;
+    cout << endl << "Takeoff report for " << name << " (Liftoff at time: "
+         << fdmex->GetState()->Getsim_time() << " 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): " << fdmex->GetPropagate()->GetAltitudeASL() << " ft. / "
+         << fdmex->GetPropagate()->GetAltitudeASLmeters() << " m  | Temperature: "
+         << fdmex->GetAtmosphere()->GetTemperature() - 459.67 << " F / "
+         << RankineToCelsius(fdmex->GetAtmosphere()->GetTemperature()) << " C]" << endl;
+    cout << "  [Velocity (KCAS): " << fdmex->GetAuxiliary()->GetVcalibratedKTS() << "]" << endl;
     TakeoffReported = true;
     break;
+  case erNone:
+    break;
   }
 }
 
@@ -809,7 +838,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)