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