src/FDM/JSBSim/models/FGLGear.cpp

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Module:       FGLGear.cpp
   4  Author:       Jon S. Berndt
   5                Norman H. Princen
   6                Bertrand Coconnier
   7  Date started: 11/18/99
   8  Purpose:      Encapsulates the landing gear elements
   9  Called by:    FGAircraft
  10
  11  ------------- Copyright (C) 1999  Jon S. Berndt (jon@jsbsim.org) -------------
  12
  13  This program is free software; you can redistribute it and/or modify it under
  14  the terms of the GNU Lesser General Public License as published by the Free Software
  15  Foundation; either version 2 of the License, or (at your option) any later
  16  version.
  17
  18  This program is distributed in the hope that it will be useful, but WITHOUT
  19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  20  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  21  details.
  22
  23  You should have received a copy of the GNU Lesser General Public License along with
  24  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  25  Place - Suite 330, Boston, MA  02111-1307, USA.
  26
  27  Further information about the GNU Lesser General Public License can also be found on
  28  the world wide web at http://www.gnu.org.
  29
  30 FUNCTIONAL DESCRIPTION
  31 --------------------------------------------------------------------------------
  32
  33 HISTORY
  34 --------------------------------------------------------------------------------
  35 11/18/99   JSB   Created
  36 01/30/01   NHP   Extended gear model to properly simulate steering and braking
  37 07/08/09   BC    Modified gear model to support large angles between aircraft and ground
  38
  39 /%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  40 INCLUDES
  41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  42
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include "FGLGear.h"
  47 #include "input_output/FGPropertyManager.h"
  48 #include "models/FGGroundReactions.h"
  49 #include "math/FGTable.h"
  50
  51 using namespace std;
  52
  53 namespace JSBSim {
  54
  55 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  56 DEFINITIONS
  57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  58
  59 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  60 GLOBAL DATA
  61 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  62
  63 static const char *IdSrc = "$Id: FGLGear.cpp,v 1.92 2011/11/10 12:06:14 jberndt Exp $";
  64 static const char *IdHdr = ID_LGEAR;
  65
  66 // Body To Structural (body frame is rotated 180 deg about Y and lengths are given in
  67 // ft instead of inches)
  68 const FGMatrix33 FGLGear::Tb2s(-1./inchtoft, 0., 0., 0., 1./inchtoft, 0., 0., 0., -1./inchtoft);
  69
  70 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  71 CLASS IMPLEMENTATION
  72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  73
  74 FGLGear::FGLGear(Element* el, FGFDMExec* fdmex, int number, const struct Inputs& inputs) :
  75   FGForce(fdmex),
  76   in(inputs),
  77   GearNumber(number),
  78   SteerAngle(0.0),
  79   Castered(false),
  80   StaticFriction(false)
  81 {
  82   Element *force_table=0;
  83   Element *dampCoeff=0;
  84   Element *dampCoeffRebound=0;
  85   string force_type="";
  86
  87   kSpring = bDamp = bDampRebound = dynamicFCoeff = staticFCoeff = rollingFCoeff = maxSteerAngle = 0;
  88   sSteerType = sBrakeGroup = sSteerType = "";
  89   isRetractable = 0;
  90   eDampType = dtLinear;
  91   eDampTypeRebound = dtLinear;
  92
  93   name = el->GetAttributeValue("name");
  94   sContactType = el->GetAttributeValue("type");
  95   if (sContactType == "BOGEY") {
  96     eContactType = ctBOGEY;
  97   } else if (sContactType == "STRUCTURE") {
  98     eContactType = ctSTRUCTURE;
  99   } else {
 100     // Unknown contact point types will be treated as STRUCTURE.
 101     eContactType = ctSTRUCTURE;
 102   }
 103
 104   // Default values for structural contact points
 105   if (eContactType == ctSTRUCTURE) {
 106     kSpring = in.EmptyWeight;
 107     bDamp = kSpring;
 108     bDampRebound = kSpring * 10;
 109     staticFCoeff = 1.0;
 110     dynamicFCoeff = 1.0;
 111   }
 112
 113   if (el->FindElement("spring_coeff"))
 114     kSpring = el->FindElementValueAsNumberConvertTo("spring_coeff", "LBS/FT");
 115   if (el->FindElement("damping_coeff")) {
 116     dampCoeff = el->FindElement("damping_coeff");
 117     if (dampCoeff->GetAttributeValue("type") == "SQUARE") {
 118       eDampType = dtSquare;
 119       bDamp   = el->FindElementValueAsNumberConvertTo("damping_coeff", "LBS/FT2/SEC2");
 120     } else {
 121       bDamp   = el->FindElementValueAsNumberConvertTo("damping_coeff", "LBS/FT/SEC");
 122     }
 123   }
 124
 125   if (el->FindElement("damping_coeff_rebound")) {
 126     dampCoeffRebound = el->FindElement("damping_coeff_rebound");
 127     if (dampCoeffRebound->GetAttributeValue("type") == "SQUARE") {
 128       eDampTypeRebound = dtSquare;
 129       bDampRebound   = el->FindElementValueAsNumberConvertTo("damping_coeff_rebound", "LBS/FT2/SEC2");
 130     } else {
 131       bDampRebound   = el->FindElementValueAsNumberConvertTo("damping_coeff_rebound", "LBS/FT/SEC");
 132     }
 133   } else {
 134     bDampRebound   = bDamp;
 135     eDampTypeRebound = eDampType;
 136   }
 137
 138   if (el->FindElement("dynamic_friction"))
 139     dynamicFCoeff = el->FindElementValueAsNumber("dynamic_friction");
 140   if (el->FindElement("static_friction"))
 141     staticFCoeff = el->FindElementValueAsNumber("static_friction");
 142   if (el->FindElement("rolling_friction"))
 143     rollingFCoeff = el->FindElementValueAsNumber("rolling_friction");
 144   if (el->FindElement("max_steer"))
 145     maxSteerAngle = el->FindElementValueAsNumberConvertTo("max_steer", "DEG");
 146   if (el->FindElement("retractable"))
 147     isRetractable = ((unsigned int)el->FindElementValueAsNumber("retractable"))>0.0?true:false;
 148
 149   GroundReactions = fdmex->GetGroundReactions();
 150   PropertyManager = fdmex->GetPropertyManager();
 151
 152   ForceY_Table = 0;
 153   force_table = el->FindElement("table");
 154   while (force_table) {
 155     force_type = force_table->GetAttributeValue("type");
 156     if (force_type == "CORNERING_COEFF") {
 157       ForceY_Table = new FGTable(PropertyManager, force_table);
 158     } else {
 159       cerr << "Undefined force table for " << name << " contact point" << endl;
 160     }
 161     force_table = el->FindNextElement("table");
 162   }
 163
 164   sBrakeGroup = el->FindElementValue("brake_group");
 165
 166   if (maxSteerAngle == 360) sSteerType = "CASTERED";
 167   else if (maxSteerAngle == 0.0) sSteerType = "FIXED";
 168   else sSteerType = "STEERABLE";
 169
 170   Element* element = el->FindElement("location");
 171   if (element) vXYZn = element->FindElementTripletConvertTo("IN");
 172   else {cerr << "No location given for contact " << name << endl; exit(-1);}
 173   SetTransformType(FGForce::tCustom);
 174
 175   element = el->FindElement("orientation");
 176   if (element && (eContactType == ctBOGEY)) {
 177     vGearOrient = element->FindElementTripletConvertTo("RAD");
 178
 179     double cp,sp,cr,sr,cy,sy;
 180
 181     cp=cos(vGearOrient(ePitch)); sp=sin(vGearOrient(ePitch));
 182     cr=cos(vGearOrient(eRoll));  sr=sin(vGearOrient(eRoll));
 183     cy=cos(vGearOrient(eYaw));   sy=sin(vGearOrient(eYaw));
 184
 185     mTGear(1,1) =  cp*cy;
 186     mTGear(2,1) =  cp*sy;
 187     mTGear(3,1) = -sp;
 188
 189     mTGear(1,2) = sr*sp*cy - cr*sy;
 190     mTGear(2,2) = sr*sp*sy + cr*cy;
 191     mTGear(3,2) = sr*cp;
 192
 193     mTGear(1,3) = cr*sp*cy + sr*sy;
 194     mTGear(2,3) = cr*sp*sy - sr*cy;
 195     mTGear(3,3) = cr*cp;
 196   }
 197   else {
 198     mTGear(1,1) = 1.;
 199     mTGear(2,2) = 1.;
 200     mTGear(3,3) = 1.;
 201   }
 202
 203   if      (sBrakeGroup == "LEFT"  ) eBrakeGrp = bgLeft;
 204   else if (sBrakeGroup == "RIGHT" ) eBrakeGrp = bgRight;
 205   else if (sBrakeGroup == "CENTER") eBrakeGrp = bgCenter;
 206   else if (sBrakeGroup == "NOSE"  ) eBrakeGrp = bgNose;
 207   else if (sBrakeGroup == "TAIL"  ) eBrakeGrp = bgTail;
 208   else if (sBrakeGroup == "NONE"  ) eBrakeGrp = bgNone;
 209   else if (sBrakeGroup.empty()    ) {eBrakeGrp = bgNone;
 210                                      sBrakeGroup = "NONE (defaulted)";}
 211   else {
 212     cerr << "Improper braking group specification in config file: "
 213          << sBrakeGroup << " is undefined." << endl;
 214   }
 215
 216   if      (sSteerType == "STEERABLE") eSteerType = stSteer;
 217   else if (sSteerType == "FIXED"    ) eSteerType = stFixed;
 218   else if (sSteerType == "CASTERED" ) {eSteerType = stCaster; Castered = true;}
 219   else if (sSteerType.empty()       ) {eSteerType = stFixed;
 220                                        sSteerType = "FIXED (defaulted)";}
 221   else {
 222     cerr << "Improper steering type specification in config file: "
 223          << sSteerType << " is undefined." << endl;
 224   }
 225
 226   GearUp = false;
 227   GearDown = true;
 228   GearPos  = 1.0;
 229   useFCSGearPos = false;
 230   Servicable = true;
 231
 232 // Add some AI here to determine if gear is located properly according to its
 233 // brake group type ??
 234
 235   WOW = lastWOW = false;
 236   ReportEnable = true;
 237   FirstContact = false;
 238   StartedGroundRun = false;
 239   TakeoffReported = LandingReported = false;
 240   LandingDistanceTraveled = TakeoffDistanceTraveled = TakeoffDistanceTraveled50ft = 0.0;
 241   MaximumStrutForce = MaximumStrutTravel = 0.0;
 242   SinkRate = GroundSpeed = 0.0;
 243
 244   vWhlVelVec.InitMatrix();
 245
 246   compressLength  = 0.0;
 247   compressSpeed   = 0.0;
 248   brakePct        = 0.0;
 249   maxCompLen      = 0.0;
 250
 251   WheelSlip = 0.0;
 252   TirePressureNorm = 1.0;
 253
 254   // Set Pacejka terms
 255
 256   Stiffness = 0.06;
 257   Shape = 2.8;
 258   Peak = staticFCoeff;
 259   Curvature = 1.03;
 260
 261   // Initialize Lagrange multipliers
 262   memset(LMultiplier, 0, sizeof(LMultiplier));
 263
 264   Debug(0);
 265 }
 266
 267 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 268
 269 FGLGear::~FGLGear()
 270 {
 271   delete ForceY_Table;
 272   Debug(1);
 273 }
 274
 275 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 276
 277 const FGColumnVector3& FGLGear::GetBodyForces(void)
 278 {
 279   double t = fdmex->GetSimTime();
 280
 281   vFn.InitMatrix();
 282
 283   if (isRetractable) ComputeRetractionState();
 284
 285   if (GearDown) {
 286     FGColumnVector3 terrainVel, dummy;
 287
 288     vLocalGear = in.Tb2l * in.vWhlBodyVec[GearNumber]; // Get local frame wheel location
 289
 290     gearLoc = in.Location.LocalToLocation(vLocalGear);
 291     // Compute the height of the theoretical location of the wheel (if strut is
 292     // not compressed) with respect to the ground level
 293     double height = gearLoc.GetContactPoint(t, contact, normal, terrainVel, dummy);
 294     vGroundNormal = in.Tec2b * normal;
 295
 296     // The height returned above is the AGL and is expressed in the Z direction
 297     // of the ECEF coordinate frame. We now need to transform this height in
 298     // actual compression of the strut (BOGEY) of in the normal direction to the
 299     // ground (STRUCTURE)
 300     double normalZ = (in.Tec2l*normal)(eZ);
 301     double LGearProj = -(mTGear.Transposed() * vGroundNormal)(eZ);
 302
 303     switch (eContactType) {
 304     case ctBOGEY:
 305       compressLength = LGearProj > 0.0 ? height * normalZ / LGearProj : 0.0;
 306       break;
 307     case ctSTRUCTURE:
 308       compressLength = height * normalZ / DotProduct(normal, normal);
 309       break;
 310     }
 311
 312     if (compressLength > 0.00) {
 313       WOW = true;
 314
 315       // The following equations use the vector to the tire contact patch
 316       // including the strut compression.
 317       FGColumnVector3 vWhlDisplVec;
 318
 319       switch(eContactType) {
 320       case ctBOGEY:
 321         vWhlDisplVec = mTGear * FGColumnVector3(0., 0., -compressLength);
 322         break;
 323       case ctSTRUCTURE:
 324         vWhlDisplVec = compressLength * vGroundNormal;
 325         break;
 326       }
 327
 328       FGColumnVector3 vWhlContactVec = in.vWhlBodyVec[GearNumber] + vWhlDisplVec;
 329       vActingXYZn = vXYZn + Tb2s * vWhlDisplVec;
 330       FGColumnVector3 vBodyWhlVel = in.PQR * vWhlContactVec;
 331       vBodyWhlVel += in.UVW - in.Tec2b * terrainVel;
 332
 333       vWhlVelVec = mTGear.Transposed() * vBodyWhlVel;
 334
 335       InitializeReporting();
 336       ComputeSteeringAngle();
 337       ComputeGroundCoordSys();
 338
 339       vLocalWhlVel = Transform().Transposed() * vBodyWhlVel;
 340
 341       if (fdmex->GetTrimStatus())
 342         compressSpeed = 0.0; // Steady state is sought during trimming
 343       else {
 344         compressSpeed = -vLocalWhlVel(eX);
 345         if (eContactType == ctBOGEY) compressSpeed /= LGearProj;
 346       }
 347
 348       ComputeVerticalStrutForce();
 349
 350       // Compute the friction coefficients in the wheel ground plane.
 351       if (eContactType == ctBOGEY) {
 352         ComputeSlipAngle();
 353         ComputeBrakeForceCoefficient();
 354         ComputeSideForceCoefficient();
 355       }
 356
 357       // Prepare the Jacobians and the Lagrange multipliers for later friction
 358       // forces calculations.
 359       ComputeJacobian(vWhlContactVec);
 360
 361     } else { // Gear is NOT compressed
 362
 363       WOW = false;
 364       compressLength = 0.0;
 365       compressSpeed = 0.0;
 366       WheelSlip = 0.0;
 367       StrutForce = 0.0;
 368
 369       // Let wheel spin down slowly
 370       vWhlVelVec(eX) -= 13.0 * in.TotalDeltaT;
 371       if (vWhlVelVec(eX) < 0.0) vWhlVelVec(eX) = 0.0;
 372
 373       // Return to neutral position between 1.0 and 0.8 gear pos.
 374       SteerAngle *= max(GetGearUnitPos()-0.8, 0.0)/0.2;
 375
 376       ResetReporting();
 377     }
 378   }
 379
 380   if (!fdmex->GetTrimStatus()) {
 381     ReportTakeoffOrLanding();
 382
 383     // Require both WOW and LastWOW to be true before checking crash conditions
 384     // to allow the WOW flag to be used in terminating a scripted run.
 385     if (WOW && lastWOW) CrashDetect();
 386
 387     lastWOW = WOW;
 388   }
 389
 390   return FGForce::GetBodyForces();
 391 }
 392
 393 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 394 // Build a local "ground" coordinate system defined by
 395 //  eX : normal to the ground
 396 //  eY : projection of the rolling direction on the ground
 397 //  eZ : projection of the sliping direction on the ground
 398
 399 void FGLGear::ComputeGroundCoordSys(void)
 400 {
 401   // Euler angles are built up to create a local frame to describe the forces
 402   // applied to the gear by the ground. Here pitch, yaw and roll do not have
 403   // any physical meaning. It is just a convenient notation.
 404   // First, "pitch" and "yaw" are determined in order to align eX with the
 405   // ground normal.
 406   if (vGroundNormal(eZ) < -1.0)
 407     vOrient(ePitch) = 0.5*M_PI;
 408   else if (1.0 < vGroundNormal(eZ))
 409     vOrient(ePitch) = -0.5*M_PI;
 410   else
 411     vOrient(ePitch) = asin(-vGroundNormal(eZ));
 412
 413   if (fabs(vOrient(ePitch)) == 0.5*M_PI)
 414     vOrient(eYaw) = 0.;
 415   else
 416     vOrient(eYaw) = atan2(vGroundNormal(eY), vGroundNormal(eX));
 417
 418   vOrient(eRoll) = 0.;
 419   UpdateCustomTransformMatrix();
 420
 421   if (eContactType == ctBOGEY) {
 422     // In the case of a bogey, the third angle "roll" is used to align the axis eY and eZ
 423     // to the rolling and sliping direction respectively.
 424     FGColumnVector3 updatedRollingAxis = Transform().Transposed() * mTGear
 425                                        * FGColumnVector3(-sin(SteerAngle), cos(SteerAngle), 0.);
 426
 427     vOrient(eRoll) = atan2(updatedRollingAxis(eY), -updatedRollingAxis(eZ));
 428     UpdateCustomTransformMatrix();
 429   }
 430 }
 431
 432 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 433
 434 void FGLGear::ComputeRetractionState(void)
 435 {
 436   double gearPos = GetGearUnitPos();
 437   if (gearPos < 0.01) {
 438     GearUp   = true;
 439     WOW      = false;
 440     GearDown = false;
 441     vWhlVelVec.InitMatrix();
 442   } else if (gearPos > 0.99) {
 443     GearDown = true;
 444     GearUp   = false;
 445   } else {
 446     GearUp   = false;
 447     GearDown = false;
 448   }
 449 }
 450
 451 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 452 // Calculate tire slip angle.
 453
 454 void FGLGear::ComputeSlipAngle(void)
 455 {
 456 // Check that the speed is non-null otherwise use the current angle
 457   if (vLocalWhlVel.Magnitude(eY,eZ) > 1E-3)
 458     WheelSlip = -atan2(vLocalWhlVel(eZ), fabs(vLocalWhlVel(eY)))*radtodeg;
 459 }
 460
 461 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 462 // Compute the steering angle in any case.
 463 // This will also make sure that animations will look right.
 464
 465 void FGLGear::ComputeSteeringAngle(void)
 466 {
 467   switch (eSteerType) {
 468   case stSteer:
 469     SteerAngle = degtorad * in.SteerPosDeg[GearNumber];
 470     break;
 471   case stFixed:
 472     SteerAngle = 0.0;
 473     break;
 474   case stCaster:
 475     if (!Castered)
 476       SteerAngle = degtorad * in.SteerPosDeg[GearNumber];
 477     else {
 478       // Check that the speed is non-null otherwise use the current angle
 479       if (vWhlVelVec.Magnitude(eX,eY) > 0.1)
 480         SteerAngle = atan2(vWhlVelVec(eY), fabs(vWhlVelVec(eX)));
 481     }
 482     break;
 483   default:
 484     cerr << "Improper steering type membership detected for this gear." << endl;
 485     break;
 486   }
 487 }
 488
 489 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 490 // Reset reporting functionality after takeoff
 491
 492 void FGLGear::ResetReporting(void)
 493 {
 494   if (in.DistanceAGL > 200.0) {
 495     FirstContact = false;
 496     StartedGroundRun = false;
 497     LandingReported = false;
 498     TakeoffReported = true;
 499     LandingDistanceTraveled = 0.0;
 500     MaximumStrutForce = MaximumStrutTravel = 0.0;
 501   }
 502 }
 503
 504 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 505
 506 void FGLGear::InitializeReporting(void)
 507 {
 508   // If this is the first time the wheel has made contact, remember some values
 509   // for later printout.
 510
 511   if (!FirstContact) {
 512     FirstContact  = true;
 513     SinkRate      =  compressSpeed;
 514     GroundSpeed   =  in.Vground;
 515     TakeoffReported = false;
 516   }
 517
 518   // If the takeoff run is starting, initialize.
 519
 520   if ((in.Vground > 0.1) &&
 521       (in.BrakePos[bgLeft] == 0) &&
 522       (in.BrakePos[bgRight] == 0) &&
 523       (in.TakeoffThrottle && !StartedGroundRun))
 524   {
 525     TakeoffDistanceTraveled = 0;
 526     TakeoffDistanceTraveled50ft = 0;
 527     StartedGroundRun = true;
 528   }
 529 }
 530
 531 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 532 // Takeoff and landing reporting functionality
 533
 534 void FGLGear::ReportTakeoffOrLanding(void)
 535 {
 536   if (FirstContact)
 537     LandingDistanceTraveled += in.Vground * in.TotalDeltaT;
 538
 539   if (StartedGroundRun) {
 540     TakeoffDistanceTraveled50ft += in.Vground * in.TotalDeltaT;
 541     if (WOW) TakeoffDistanceTraveled += in.Vground * in.TotalDeltaT;
 542   }
 543
 544   if ( ReportEnable
 545        && in.Vground <= 0.05
 546        && !LandingReported
 547        && in.WOW)
 548   {
 549     if (debug_lvl > 0) Report(erLand);
 550   }
 551
 552   if ( ReportEnable
 553        && !TakeoffReported
 554        && (in.DistanceAGL - vLocalGear(eZ)) > 50.0
 555        && !in.WOW)
 556   {
 557     if (debug_lvl > 0) Report(erTakeoff);
 558   }
 559
 560   if (lastWOW != WOW) PutMessage("GEAR_CONTACT: " + name, WOW);
 561 }
 562
 563 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 564 // Crash detection logic (really out-of-bounds detection)
 565
 566 void FGLGear::CrashDetect(void)
 567 {
 568   if ( (compressLength > 500.0 ||
 569       vFn.Magnitude() > 100000000.0 ||
 570       GetMoments().Magnitude() > 5000000000.0 ||
 571       SinkRate > 1.4666*30 ) && !fdmex->IntegrationSuspended())
 572   {
 573     PutMessage("Crash Detected: Simulation FREEZE.");
 574     fdmex->SuspendIntegration();
 575   }
 576 }
 577
 578 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 579 // The following needs work regarding friction coefficients and braking and
 580 // steering The BrakeFCoeff formula assumes that an anti-skid system is used.
 581 // It also assumes that we won't be turning and braking at the same time.
 582 // Will fix this later.
 583 // [JSB] The braking force coefficients include normal rolling coefficient +
 584 // a percentage of the static friction coefficient based on braking applied.
 585
 586 void FGLGear::ComputeBrakeForceCoefficient(void)
 587 {
 588   switch (eBrakeGrp) {
 589   case bgLeft:
 590     BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgLeft]) +
 591                      staticFCoeff * in.BrakePos[bgLeft] );
 592     break;
 593   case bgRight:
 594     BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgRight]) +
 595                      staticFCoeff * in.BrakePos[bgRight] );
 596     break;
 597   case bgCenter:
 598     BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgCenter]) +
 599                      staticFCoeff * in.BrakePos[bgCenter] );
 600     break;
 601   case bgNose:
 602     BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgCenter]) +
 603                      staticFCoeff * in.BrakePos[bgCenter] );
 604     break;
 605   case bgTail:
 606     BrakeFCoeff =  ( rollingFCoeff * (1.0 - in.BrakePos[bgCenter]) +
 607                      staticFCoeff * in.BrakePos[bgCenter] );
 608     break;
 609   case bgNone:
 610     BrakeFCoeff =  rollingFCoeff;
 611     break;
 612   default:
 613     cerr << "Improper brake group membership detected for this gear." << endl;
 614     break;
 615   }
 616 }
 617
 618 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 619 // Compute the sideforce coefficients using Pacejka's Magic Formula.
 620 //
 621 //   y(x) = D sin {C arctan [Bx - E(Bx - arctan Bx)]}
 622 //
 623 // Where: B = Stiffness Factor (0.06, here)
 624 //        C = Shape Factor (2.8, here)
 625 //        D = Peak Factor (0.8, here)
 626 //        E = Curvature Factor (1.03, here)
 627
 628 void FGLGear::ComputeSideForceCoefficient(void)
 629 {
 630   if (ForceY_Table) {
 631     FCoeff = ForceY_Table->GetValue(WheelSlip);
 632   } else {
 633     double StiffSlip = Stiffness*WheelSlip;
 634     FCoeff = Peak * sin(Shape*atan(StiffSlip - Curvature*(StiffSlip - atan(StiffSlip))));
 635   }
 636 }
 637
 638 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 639 // Compute the vertical force on the wheel using square-law damping (per comment
 640 // in paper AIAA-2000-4303 - see header prologue comments). We might consider
 641 // allowing for both square and linear damping force calculation. Also need to
 642 // possibly give a "rebound damping factor" that differs from the compression
 643 // case.
 644
 645 void FGLGear::ComputeVerticalStrutForce(void)
 646 {
 647   double springForce = 0;
 648   double dampForce = 0;
 649
 650   springForce = -compressLength * kSpring;
 651
 652   if (compressSpeed >= 0.0) {
 653
 654     if (eDampType == dtLinear)   dampForce = -compressSpeed * bDamp;
 655     else         dampForce = -compressSpeed * compressSpeed * bDamp;
 656
 657   } else {
 658
 659     if (eDampTypeRebound == dtLinear)
 660       dampForce   = -compressSpeed * bDampRebound;
 661     else
 662       dampForce   =  compressSpeed * compressSpeed * bDampRebound;
 663
 664   }
 665
 666   StrutForce = min(springForce + dampForce, (double)0.0);
 667
 668   // The reaction force of the wheel is always normal to the ground
 669   switch (eContactType) {
 670   case ctBOGEY:
 671     // Project back the strut force in the local coordinate frame of the ground
 672     vFn(eX) = StrutForce / (mTGear.Transposed()*vGroundNormal)(eZ);
 673     break;
 674   case ctSTRUCTURE:
 675     vFn(eX) = -StrutForce;
 676     break;
 677   }
 678
 679   // Remember these values for reporting
 680   MaximumStrutForce = max(MaximumStrutForce, fabs(StrutForce));
 681   MaximumStrutTravel = max(MaximumStrutTravel, fabs(compressLength));
 682 }
 683
 684 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 685
 686 double FGLGear::GetGearUnitPos(void)
 687 {
 688   // hack to provide backward compatibility to gear/gear-pos-norm property
 689   if( useFCSGearPos || in.FCSGearPos != 1.0 ) {
 690     useFCSGearPos = true;
 691     return in.FCSGearPos;
 692   }
 693   return GearPos;
 694 }
 695
 696 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 697 // Compute the jacobian entries for the friction forces resolution later
 698 // in FGPropagate
 699
 700 void FGLGear::ComputeJacobian(const FGColumnVector3& vWhlContactVec)
 701 {
 702   // When the point of contact is moving, dynamic friction is used
 703   // This type of friction is limited to ctSTRUCTURE elements because their
 704   // friction coefficient is the same in every directions
 705   if ((eContactType == ctSTRUCTURE) && (vLocalWhlVel.Magnitude(eY,eZ) > 1E-3)) {
 706     FGColumnVector3 velocityDirection = vLocalWhlVel;
 707
 708     StaticFriction = false;
 709
 710     velocityDirection(eX) = 0.;
 711     velocityDirection.Normalize();
 712
 713     LMultiplier[ftDynamic].ForceJacobian = Transform()*velocityDirection;
 714     LMultiplier[ftDynamic].MomentJacobian = vWhlContactVec * LMultiplier[ftDynamic].ForceJacobian;
 715     LMultiplier[ftDynamic].Max = 0.;
 716     LMultiplier[ftDynamic].Min = -fabs(dynamicFCoeff * vFn(eX));
 717
 718     // The Lagrange multiplier value obtained from the previous iteration is kept
 719     // This is supposed to accelerate the convergence of the projected Gauss-Seidel
 720     // algorithm. The code just below is to make sure that the initial value
 721     // is consistent with the current friction coefficient and normal reaction.
 722     LMultiplier[ftDynamic].value = Constrain(LMultiplier[ftDynamic].Min, LMultiplier[ftDynamic].value, LMultiplier[ftDynamic].Max);
 723
 724     GroundReactions->RegisterLagrangeMultiplier(&LMultiplier[ftDynamic]);
 725   }
 726   else {
 727     // Static friction is used for ctSTRUCTURE when the contact point is not moving.
 728     // It is always used for ctBOGEY elements because the friction coefficients
 729     // of a tyre depend on the direction of the movement (roll & side directions).
 730     // This cannot be handled properly by the so-called "dynamic friction".
 731     StaticFriction = true;
 732
 733     LMultiplier[ftRoll].ForceJacobian = Transform()*FGColumnVector3(0.,1.,0.);
 734     LMultiplier[ftSide].ForceJacobian = Transform()*FGColumnVector3(0.,0.,1.);
 735     LMultiplier[ftRoll].MomentJacobian = vWhlContactVec * LMultiplier[ftRoll].ForceJacobian;
 736     LMultiplier[ftSide].MomentJacobian = vWhlContactVec * LMultiplier[ftSide].ForceJacobian;
 737
 738     switch(eContactType) {
 739     case ctBOGEY:
 740       LMultiplier[ftRoll].Max = fabs(BrakeFCoeff * vFn(eX));
 741       LMultiplier[ftSide].Max = fabs(FCoeff * vFn(eX));
 742       break;
 743     case ctSTRUCTURE:
 744       LMultiplier[ftRoll].Max = fabs(staticFCoeff * vFn(eX));
 745       LMultiplier[ftSide].Max = fabs(staticFCoeff * vFn(eX));
 746       break;
 747     }
 748
 749     LMultiplier[ftRoll].Min = -LMultiplier[ftRoll].Max;
 750     LMultiplier[ftSide].Min = -LMultiplier[ftSide].Max;
 751
 752     // The Lagrange multiplier value obtained from the previous iteration is kept
 753     // This is supposed to accelerate the convergence of the projected Gauss-Seidel
 754     // algorithm. The code just below is to make sure that the initial value
 755     // is consistent with the current friction coefficient and normal reaction.
 756     LMultiplier[ftRoll].value = Constrain(LMultiplier[ftRoll].Min, LMultiplier[ftRoll].value, LMultiplier[ftRoll].Max);
 757     LMultiplier[ftSide].value = Constrain(LMultiplier[ftSide].Min, LMultiplier[ftSide].value, LMultiplier[ftSide].Max);
 758
 759     GroundReactions->RegisterLagrangeMultiplier(&LMultiplier[ftRoll]);
 760     GroundReactions->RegisterLagrangeMultiplier(&LMultiplier[ftSide]);
 761   }
 762 }
 763
 764 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 765 // This routine is called after the Lagrange multiplier has been computed in
 766 // the FGAccelerations class. The friction forces of the landing gear are then
 767 // updated accordingly.
 768 void FGLGear::UpdateForces(void)
 769 {
 770   if (StaticFriction) {
 771     vFn(eY) = LMultiplier[ftRoll].value;
 772     vFn(eZ) = LMultiplier[ftSide].value;
 773   }
 774   else
 775     vFn += LMultiplier[ftDynamic].value * (Transform ().Transposed() * LMultiplier[ftDynamic].ForceJacobian);
 776 }
 777
 778 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 779
 780 void FGLGear::bind(void)
 781 {
 782   string property_name;
 783   string base_property_name;
 784
 785   switch(eContactType) {
 786   case ctBOGEY:
 787     base_property_name = CreateIndexedPropertyName("gear/unit", GearNumber);
 788     break;
 789   case ctSTRUCTURE:
 790     base_property_name = CreateIndexedPropertyName("contact/unit", GearNumber);
 791     break;
 792   default:
 793     return;
 794   }
 795
 796   property_name = base_property_name + "/WOW";
 797   PropertyManager->Tie( property_name.c_str(), &WOW );
 798   property_name = base_property_name + "/z-position";
 799   PropertyManager->Tie( property_name.c_str(), (FGForce*)this,
 800                           &FGForce::GetLocationZ, &FGForce::SetLocationZ);
 801   property_name = base_property_name + "/compression-ft";
 802   PropertyManager->Tie( property_name.c_str(), &compressLength );
 803   property_name = base_property_name + "/static_friction_coeff";
 804   PropertyManager->Tie( property_name.c_str(), &staticFCoeff );
 805   property_name = base_property_name + "/dynamic_friction_coeff";
 806   PropertyManager->Tie( property_name.c_str(), &dynamicFCoeff );
 807
 808   if (eContactType == ctBOGEY) {
 809     property_name = base_property_name + "/slip-angle-deg";
 810     PropertyManager->Tie( property_name.c_str(), &WheelSlip );
 811     property_name = base_property_name + "/wheel-speed-fps";
 812     PropertyManager->Tie( property_name.c_str(), (FGLGear*)this,
 813                           &FGLGear::GetWheelRollVel);
 814     property_name = base_property_name + "/side_friction_coeff";
 815     PropertyManager->Tie( property_name.c_str(), &FCoeff );
 816     property_name = base_property_name + "/rolling_friction_coeff";
 817     PropertyManager->Tie( property_name.c_str(), &rollingFCoeff );
 818
 819     if (eSteerType == stCaster) {
 820       property_name = base_property_name + "/steering-angle-deg";
 821       PropertyManager->Tie( property_name.c_str(), this, &FGLGear::GetSteerAngleDeg );
 822       property_name = base_property_name + "/castered";
 823       PropertyManager->Tie( property_name.c_str(), &Castered);
 824     }
 825   }
 826
 827   if( isRetractable ) {
 828     property_name = base_property_name + "/pos-norm";
 829     PropertyManager->Tie( property_name.c_str(), &GearPos );
 830   }
 831 }
 832
 833 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 834
 835 void FGLGear::Report(ReportType repType)
 836 {
 837   if (fabs(TakeoffDistanceTraveled) < 0.001) return; // Don't print superfluous reports
 838
 839   switch(repType) {
 840   case erLand:
 841     cout << endl << "Touchdown report for " << name << " (WOW at time: "
 842          << fdmex->GetSimTime() << " seconds)" << endl;
 843     cout << "  Sink rate at contact:  " << SinkRate                << " fps,    "
 844                                 << SinkRate*0.3048          << " mps"     << endl;
 845     cout << "  Contact ground speed:  " << GroundSpeed*.5925       << " knots,  "
 846                                 << GroundSpeed*0.3048       << " mps"     << endl;
 847     cout << "  Maximum contact force: " << MaximumStrutForce       << " lbs,    "
 848                                 << MaximumStrutForce*4.448  << " Newtons" << endl;
 849     cout << "  Maximum strut travel:  " << MaximumStrutTravel*12.0 << " inches, "
 850                                 << MaximumStrutTravel*30.48 << " cm"      << endl;
 851     cout << "  Distance traveled:     " << LandingDistanceTraveled        << " ft,     "
 852                                 << LandingDistanceTraveled*0.3048  << " meters"  << endl;
 853     LandingReported = true;
 854     break;
 855   case erTakeoff:
 856     cout << endl << "Takeoff report for " << name << " (Liftoff at time: "
 857         << fdmex->GetSimTime() << " seconds)" << endl;
 858     cout << "  Distance traveled:                " << TakeoffDistanceTraveled
 859          << " ft,     " << TakeoffDistanceTraveled*0.3048  << " meters"  << endl;
 860     cout << "  Distance traveled (over 50'):     " << TakeoffDistanceTraveled50ft
 861          << " ft,     " << TakeoffDistanceTraveled50ft*0.3048 << " meters" << endl;
 862     cout << "  [Altitude (ASL): " << in.DistanceASL << " ft. / "
 863          << in.DistanceASL*FGJSBBase::fttom << " m  | Temperature: "
 864          << in.Temperature - 459.67 << " F / "
 865          << RankineToCelsius(in.Temperature) << " C]" << endl;
 866     cout << "  [Velocity (KCAS): " << in.VcalibratedKts << "]" << endl;
 867     TakeoffReported = true;
 868     break;
 869   case erNone:
 870     break;
 871   }
 872 }
 873
 874 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 875 //    The bitmasked value choices are as follows:
 876 //    unset: In this case (the default) JSBSim would only print
 877 //       out the normally expected messages, essentially echoing
 878 //       the config files as they are read. If the environment
 879 //       variable is not set, debug_lvl is set to 1 internally
 880 //    0: This requests JSBSim not to output any messages
 881 //       whatsoever.
 882 //    1: This value explicity requests the normal JSBSim
 883 //       startup messages
 884 //    2: This value asks for a message to be printed out when
 885 //       a class is instantiated
 886 //    4: When this value is set, a message is displayed when a
 887 //       FGModel object executes its Run() method
 888 //    8: When this value is set, various runtime state variables
 889 //       are printed out periodically
 890 //    16: When set various parameters are sanity checked and
 891 //       a message is printed out when they go out of bounds
 892
 893 void FGLGear::Debug(int from)
 894 {
 895   if (debug_lvl <= 0) return;
 896
 897   if (debug_lvl & 1) { // Standard console startup message output
 898     if (from == 0) { // Constructor - loading and initialization
 899       cout << "    " << sContactType << " " << name          << endl;
 900       cout << "      Location: "         << vXYZn          << endl;
 901       cout << "      Spring Constant:  " << kSpring       << endl;
 902
 903       if (eDampType == dtLinear)
 904         cout << "      Damping Constant: " << bDamp << " (linear)" << endl;
 905       else
 906         cout << "      Damping Constant: " << bDamp << " (square law)" << endl;
 907
 908       if (eDampTypeRebound == dtLinear)
 909         cout << "      Rebound Damping Constant: " << bDampRebound << " (linear)" << endl;
 910       else
 911         cout << "      Rebound Damping Constant: " << bDampRebound << " (square law)" << endl;
 912
 913       cout << "      Dynamic Friction: " << dynamicFCoeff << endl;
 914       cout << "      Static Friction:  " << staticFCoeff  << endl;
 915       if (eContactType == ctBOGEY) {
 916         cout << "      Rolling Friction: " << rollingFCoeff << endl;
 917         cout << "      Steering Type:    " << sSteerType    << endl;
 918         cout << "      Grouping:         " << sBrakeGroup   << endl;
 919         cout << "      Max Steer Angle:  " << maxSteerAngle << endl;
 920         cout << "      Retractable:      " << isRetractable  << endl;
 921       }
 922     }
 923   }
 924   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 925     if (from == 0) cout << "Instantiated: FGLGear" << endl;
 926     if (from == 1) cout << "Destroyed:    FGLGear" << endl;
 927   }
 928   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 929   }
 930   if (debug_lvl & 8 ) { // Runtime state variables
 931   }
 932   if (debug_lvl & 16) { // Sanity checking
 933   }
 934   if (debug_lvl & 64) {
 935     if (from == 0) { // Constructor
 936       cout << IdSrc << endl;
 937       cout << IdHdr << endl;
 938     }
 939   }
 940 }
 941
 942 } // namespace JSBSim
 943