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