src/FDM/JSBSim/FGTrimAxis.cpp

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Header:       FGTrimAxis.cpp
   4  Author:       Tony Peden
   5  Date started: 7/3/00
   6
   7  --------- Copyright (C) 1999  Anthony K. Peden (apeden@earthlink.net) ---------
   8
   9  This program is free software; you can redistribute it and/or modify it under
  10  the terms of the GNU General Public License as published by the Free Software
  11  Foundation; either version 2 of the License, or (at your option) any later
  12  version.
  13
  14  This program is distributed in the hope that it will be useful, but WITHOUT
  15  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  16  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  17  details.
  18
  19  You should have received a copy of the GNU General Public License along with
  20  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  21  Place - Suite 330, Boston, MA  02111-1307, USA.
  22
  23  Further information about the GNU General Public License can also be found on
  24  the world wide web at http://www.gnu.org.
  25
  26
  27  HISTORY
  28 --------------------------------------------------------------------------------
  29 7/3/00   TP   Created
  30
  31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  32 INCLUDES
  33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  34
  35 #ifdef _MSC_VER
  36 #  pragma warning (disable : 4786)
  37 #endif
  38
  39 #include <string>
  40 #include <stdlib.h>
  41
  42 #include "FGFDMExec.h"
  43 #include "FGAtmosphere.h"
  44 #include "FGInitialCondition.h"
  45 #include "FGTrimAxis.h"
  46 #include "FGAircraft.h"
  47 #include "FGPropulsion.h"
  48 #include "FGAerodynamics.h"
  49
  50
  51 namespace JSBSim {
  52
  53 static const char *IdSrc = "$Id$";
  54 static const char *IdHdr = ID_TRIMAXIS;
  55
  56 /*****************************************************************************/
  57
  58 FGTrimAxis::FGTrimAxis(FGFDMExec* fdex, FGInitialCondition* ic, State st,
  59                        Control ctrl) {
  60
  61   fdmex=fdex;
  62   fgic=ic;
  63   state=st;
  64   control=ctrl;
  65   max_iterations=10;
  66   control_value=0;
  67   its_to_stable_value=0;
  68   total_iterations=0;
  69   total_stability_iterations=0;
  70   state_convert=1.0;
  71   control_convert=1.0;
  72   state_value=0;
  73   state_target=0;
  74   switch(state) {
  75     case tUdot: tolerance = DEFAULT_TOLERANCE; break;
  76     case tVdot: tolerance = DEFAULT_TOLERANCE; break;
  77     case tWdot: tolerance = DEFAULT_TOLERANCE; break;
  78     case tQdot: tolerance = DEFAULT_TOLERANCE / 10; break;
  79     case tPdot: tolerance = DEFAULT_TOLERANCE / 10; break;
  80     case tRdot: tolerance = DEFAULT_TOLERANCE / 10; break;
  81     case tHmgt: tolerance = 0.01; break;
  82     case  tNlf: state_target=1.0; tolerance = 1E-5; break;
  83     case tAll: break;
  84   }
  85
  86   solver_eps=tolerance;
  87   switch(control) {
  88   case tThrottle:
  89     control_min=0;
  90     control_max=1;
  91     control_value=0.5;
  92     break;
  93   case tBeta:
  94     control_min=-30*degtorad;
  95     control_max=30*degtorad;
  96     control_convert=radtodeg;
  97     break;
  98   case tAlpha:
  99     control_min=fdmex->GetAerodynamics()->GetAlphaCLMin();
 100     control_max=fdmex->GetAerodynamics()->GetAlphaCLMax();
 101     if(control_max <= control_min) {
 102       control_max=20*degtorad;
 103       control_min=-5*degtorad;
 104     }
 105     control_value= (control_min+control_max)/2;
 106     control_convert=radtodeg;
 107     solver_eps=tolerance/100;
 108     break;
 109   case tPitchTrim:
 110   case tElevator:
 111   case tRollTrim:
 112   case tAileron:
 113   case tYawTrim:
 114   case tRudder:
 115     control_min=-1;
 116     control_max=1;
 117     state_convert=radtodeg;
 118     solver_eps=tolerance/100;
 119     break;
 120   case tAltAGL:
 121     control_min=0;
 122     control_max=30;
 123     control_value=fdmex->GetPosition()->GetDistanceAGL();
 124     solver_eps=tolerance/100;
 125     break;
 126   case tTheta:
 127     control_min=fdmex->GetRotation()->Gettht() - 5*degtorad;
 128     control_max=fdmex->GetRotation()->Gettht() + 5*degtorad;
 129     state_convert=radtodeg;
 130     break;
 131   case tPhi:
 132     control_min=fdmex->GetRotation()->Getphi() - 30*degtorad;
 133     control_max=fdmex->GetRotation()->Getphi() + 30*degtorad;
 134     state_convert=radtodeg;
 135     control_convert=radtodeg;
 136     break;
 137   case tGamma:
 138     solver_eps=tolerance/100;
 139     control_min=-80*degtorad;
 140     control_max=80*degtorad;
 141     control_convert=radtodeg;
 142     break;
 143   case tHeading:
 144     control_min=fdmex->GetRotation()->Getpsi() - 30*degtorad;
 145     control_max=fdmex->GetRotation()->Getpsi() + 30*degtorad;
 146     state_convert=radtodeg;
 147     break;
 148   }
 149
 150
 151   Debug(0);
 152 }
 153
 154 /*****************************************************************************/
 155
 156 FGTrimAxis::~FGTrimAxis(void)
 157 {
 158   Debug(1);
 159 }
 160
 161 /*****************************************************************************/
 162
 163 void FGTrimAxis::getState(void) {
 164   switch(state) {
 165   case tUdot: state_value=fdmex->GetTranslation()->GetUVWdot(1)-state_target; break;
 166   case tVdot: state_value=fdmex->GetTranslation()->GetUVWdot(2)-state_target; break;
 167   case tWdot: state_value=fdmex->GetTranslation()->GetUVWdot(3)-state_target; break;
 168   case tQdot: state_value=fdmex->GetRotation()->GetPQRdot(2)-state_target;break;
 169   case tPdot: state_value=fdmex->GetRotation()->GetPQRdot(1)-state_target; break;
 170   case tRdot: state_value=fdmex->GetRotation()->GetPQRdot(3)-state_target; break;
 171   case tHmgt: state_value=computeHmgt()-state_target; break;
 172   case tNlf:  state_value=fdmex->GetAircraft()->GetNlf()-state_target; break;
 173   case tAll: break;
 174   }
 175 }
 176
 177 /*****************************************************************************/
 178
 179 //States are not settable
 180
 181 void FGTrimAxis::getControl(void) {
 182   switch(control) {
 183   case tThrottle:  control_value=fdmex->GetFCS()->GetThrottleCmd(0); break;
 184   case tBeta:      control_value=fdmex->GetTranslation()->Getalpha(); break;
 185   case tAlpha:     control_value=fdmex->GetTranslation()->Getbeta();  break;
 186   case tPitchTrim: control_value=fdmex->GetFCS() -> GetPitchTrimCmd(); break;
 187   case tElevator:  control_value=fdmex->GetFCS() -> GetDeCmd(); break;
 188   case tRollTrim:
 189   case tAileron:   control_value=fdmex->GetFCS() -> GetDaCmd(); break;
 190   case tYawTrim:
 191   case tRudder:    control_value=fdmex->GetFCS() -> GetDrCmd(); break;
 192   case tAltAGL:    control_value=fdmex->GetPosition()->GetDistanceAGL();break;
 193   case tTheta:     control_value=fdmex->GetRotation()->Gettht(); break;
 194   case tPhi:       control_value=fdmex->GetRotation()->Getphi(); break;
 195   case tGamma:     control_value=fdmex->GetPosition()->GetGamma();break;
 196   case tHeading:   control_value=fdmex->GetRotation()->Getpsi(); break;
 197   }
 198 }
 199
 200 /*****************************************************************************/
 201
 202 double FGTrimAxis::computeHmgt(void) {
 203   double diff;
 204
 205   diff   = fdmex->GetRotation()->Getpsi() -
 206              fdmex->GetPosition()->GetGroundTrack();
 207
 208   if( diff < -M_PI ) {
 209      return (diff + 2*M_PI);
 210   } else if( diff > M_PI ) {
 211      return (diff - 2*M_PI);
 212   } else {
 213      return diff;
 214   }
 215
 216 }
 217
 218 /*****************************************************************************/
 219
 220
 221 void FGTrimAxis::setControl(void) {
 222   switch(control) {
 223   case tThrottle:  setThrottlesPct(); break;
 224   case tBeta:      fgic->SetBetaRadIC(control_value); break;
 225   case tAlpha:     fgic->SetAlphaRadIC(control_value);  break;
 226   case tPitchTrim: fdmex->GetFCS()->SetPitchTrimCmd(control_value); break;
 227   case tElevator:  fdmex->GetFCS()->SetDeCmd(control_value); break;
 228   case tRollTrim:
 229   case tAileron:   fdmex->GetFCS()->SetDaCmd(control_value); break;
 230   case tYawTrim:
 231   case tRudder:    fdmex->GetFCS()->SetDrCmd(control_value); break;
 232   case tAltAGL:    fgic->SetAltitudeAGLFtIC(control_value); break;
 233   case tTheta:     fgic->SetPitchAngleRadIC(control_value); break;
 234   case tPhi:       fgic->SetRollAngleRadIC(control_value); break;
 235   case tGamma:     fgic->SetFlightPathAngleRadIC(control_value); break;
 236   case tHeading:   fgic->SetTrueHeadingRadIC(control_value); break;
 237   }
 238 }
 239
 240
 241
 242
 243
 244 /*****************************************************************************/
 245
 246 // the aircraft center of rotation is no longer the cg once the gear
 247 // contact the ground so the altitude needs to be changed when pitch
 248 // and roll angle are adjusted.  Instead of attempting to calculate the
 249 // new center of rotation, pick a gear unit as a reference and use its
 250 // location vector to calculate the new height change. i.e. new altitude =
 251 // earth z component of that vector (which is in body axes )
 252 void FGTrimAxis::SetThetaOnGround(double ff) {
 253   int center,i,ref;
 254
 255   // favor an off-center unit so that the same one can be used for both
 256   // pitch and roll.  An on-center unit is used (for pitch)if that's all
 257   // that's in contact with the ground.
 258   i=0; ref=-1; center=-1;
 259   while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
 260     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) {
 261       if(fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01)
 262         ref=i;
 263       else
 264         center=i;
 265     }
 266     i++;
 267   }
 268   if((ref < 0) && (center >= 0)) {
 269     ref=center;
 270   }
 271   cout << "SetThetaOnGround ref gear: " << ref << endl;
 272   if(ref >= 0) {
 273     double sp=fdmex->GetRotation()->GetSinphi();
 274     double cp=fdmex->GetRotation()->GetCosphi();
 275     double lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
 276     double ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
 277     double lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
 278     double hagl = -1*lx*sin(ff) +
 279                     ly*sp*cos(ff) +
 280                     lz*cp*cos(ff);
 281
 282     fgic->SetAltitudeAGLFtIC(hagl);
 283     cout << "SetThetaOnGround new alt: " << hagl << endl;
 284   }
 285   fgic->SetPitchAngleRadIC(ff);
 286   cout << "SetThetaOnGround new theta: " << ff << endl;
 287 }
 288
 289 /*****************************************************************************/
 290
 291 bool FGTrimAxis::initTheta(void) {
 292   int i,N,iAft, iForward;
 293   double zAft,zForward,zDiff,theta;
 294   bool level;
 295   double saveAlt;
 296
 297   saveAlt=fgic->GetAltitudeAGLFtIC();
 298   fgic->SetAltitudeAGLFtIC(100);
 299
 300
 301   N=fdmex->GetGroundReactions()->GetNumGearUnits();
 302
 303   //find the first wheel unit forward of the cg
 304   //the list is short so a simple linear search is fine
 305   for( i=0; i<N; i++ ) {
 306     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) > 0 ) {
 307         iForward=i;
 308         break;
 309     }
 310   }
 311   //now find the first wheel unit aft of the cg
 312   for( i=0; i<N; i++ ) {
 313     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) < 0 ) {
 314         iAft=i;
 315         break;
 316     }
 317   }
 318
 319   // now adjust theta till the wheels are the same distance from the ground
 320   zAft=fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear(3);
 321   zForward=fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear(3);
 322   zDiff = zForward - zAft;
 323   level=false;
 324   theta=fgic->GetPitchAngleDegIC();
 325   while(!level && (i < 100)) {
 326           theta+=2.0*zDiff;
 327           fgic->SetPitchAngleDegIC(theta);
 328           fdmex->RunIC();
 329           zAft=fdmex->GetGroundReactions()->GetGearUnit(iAft)->GetLocalGear(3);
 330     zForward=fdmex->GetGroundReactions()->GetGearUnit(iForward)->GetLocalGear(3);
 331     zDiff = zForward - zAft;
 332           //cout << endl << theta << "  " << zDiff << endl;
 333           //cout << "0: " << fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear() << endl;
 334           //cout << "1: " << fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear() << endl;
 335           if(fabs(zDiff ) < 0.1)
 336               level=true;
 337           i++;
 338   }
 339   //cout << i << endl;
 340   if (debug_lvl > 0) {
 341       cout << "    Initial Theta: " << fdmex->GetRotation()->Gettht()*radtodeg << endl;
 342       cout << "    Used gear unit " << iAft << " as aft and " << iForward << " as forward" << endl;
 343   }
 344   control_min=(theta+5)*degtorad;
 345   control_max=(theta-5)*degtorad;
 346   fgic->SetAltitudeAGLFtIC(saveAlt);
 347   if(i < 100)
 348     return true;
 349   else
 350     return false;
 351 }
 352
 353 /*****************************************************************************/
 354
 355 void FGTrimAxis::SetPhiOnGround(double ff) {
 356   int i,ref;
 357
 358   i=0; ref=-1;
 359   //must have an off-center unit here
 360   while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
 361     if( (fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) &&
 362       (fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01))
 363         ref=i;
 364     i++;
 365   }
 366   if(ref >= 0) {
 367     double st=fdmex->GetRotation()->GetSintht();
 368     double ct=fdmex->GetRotation()->GetCostht();
 369     double lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
 370     double ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
 371     double lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
 372     double hagl = -1*lx*st +
 373                     ly*sin(ff)*ct +
 374                     lz*cos(ff)*ct;
 375
 376     fgic->SetAltitudeAGLFtIC(hagl);
 377   }
 378   fgic->SetRollAngleRadIC(ff);
 379
 380 }
 381
 382 /*****************************************************************************/
 383
 384 void FGTrimAxis::Run(void) {
 385
 386   double last_state_value;
 387   int i;
 388   setControl();
 389   //cout << "FGTrimAxis::Run: " << control_value << endl;
 390   i=0;
 391   bool stable=false;
 392   while(!stable) {
 393     i++;
 394     last_state_value=state_value;
 395     fdmex->RunIC();
 396     getState();
 397     if(i > 1) {
 398       if((fabs(last_state_value - state_value) < tolerance) || (i >= 100) )
 399         stable=true;
 400     }
 401   }
 402
 403   its_to_stable_value=i;
 404   total_stability_iterations+=its_to_stable_value;
 405   total_iterations++;
 406 }
 407
 408 /*****************************************************************************/
 409
 410 void FGTrimAxis::setThrottlesPct(void) {
 411   double tMin,tMax;
 412   for(unsigned i=0;i<fdmex->GetPropulsion()->GetNumEngines();i++) {
 413       tMin=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMin();
 414       tMax=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMax();
 415       //cout << "setThrottlespct: " << i << ", " << control_min << ", " << control_max << ", " << control_value;
 416       fdmex->GetFCS()->SetThrottleCmd(i,tMin+control_value*(tMax-tMin));
 417       //cout << "setThrottlespct: " << fdmex->GetFCS()->GetThrottleCmd(i) << endl;
 418       fdmex->RunIC(); //apply throttle change
 419       fdmex->GetPropulsion()->GetSteadyState();
 420   }
 421 }
 422
 423 /*****************************************************************************/
 424
 425 void FGTrimAxis::AxisReport(void) {
 426
 427   char out[80];
 428   sprintf(out,"  %20s: %6.2f %5s: %9.2e Tolerance: %3.0e\n",
 429            GetControlName().c_str(), GetControl()*control_convert,
 430            GetStateName().c_str(), GetState()+state_target, GetTolerance());
 431   cout << out;
 432
 433 }
 434
 435 /*****************************************************************************/
 436
 437 double FGTrimAxis::GetAvgStability( void ) {
 438   if(total_iterations > 0) {
 439     return double(total_stability_iterations)/double(total_iterations);
 440   }
 441   return 0;
 442 }
 443
 444 /*****************************************************************************/
 445 //    The bitmasked value choices are as follows:
 446 //    unset: In this case (the default) JSBSim would only print
 447 //       out the normally expected messages, essentially echoing
 448 //       the config files as they are read. If the environment
 449 //       variable is not set, debug_lvl is set to 1 internally
 450 //    0: This requests JSBSim not to output any messages
 451 //       whatsoever.
 452 //    1: This value explicity requests the normal JSBSim
 453 //       startup messages
 454 //    2: This value asks for a message to be printed out when
 455 //       a class is instantiated
 456 //    4: When this value is set, a message is displayed when a
 457 //       FGModel object executes its Run() method
 458 //    8: When this value is set, various runtime state variables
 459 //       are printed out periodically
 460 //    16: When set various parameters are sanity checked and
 461 //       a message is printed out when they go out of bounds
 462
 463 void FGTrimAxis::Debug(int from)
 464 {
 465
 466   if (debug_lvl <= 0) return;
 467   if (debug_lvl & 1 ) { // Standard console startup message output
 468     if (from == 0) { // Constructor
 469
 470     }
 471   }
 472   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 473     if (from == 0) cout << "Instantiated: FGTrimAxis" << endl;
 474     if (from == 1) cout << "Destroyed:    FGTrimAxis" << endl;
 475   }
 476   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 477   }
 478   if (debug_lvl & 8 ) { // Runtime state variables
 479   }
 480   if (debug_lvl & 16) { // Sanity checking
 481   }
 482   if (debug_lvl & 64) {
 483     if (from == 0) { // Constructor
 484       cout << IdSrc << endl;
 485       cout << IdHdr << endl;
 486     }
 487   }
 488 }
 489 }