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