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     switch(state) {
  66     case tUdot: tolerance = DEFAULT_TOLERANCE; break;
  67     case tVdot: tolerance = DEFAULT_TOLERANCE; break;
  68     case tWdot: tolerance = DEFAULT_TOLERANCE; break;
  69     case tQdot: tolerance = DEFAULT_TOLERANCE / 10; break;
  70     case tPdot: tolerance = DEFAULT_TOLERANCE / 10; break;
  71     case tRdot: tolerance = DEFAULT_TOLERANCE / 10; break;
  72     case tHmgt: tolerance = 0.01; break;
  73   }
  74
  75   solver_eps=tolerance;
  76   switch(control) {
  77   case tThrottle:
  78     control_min=0;
  79     control_max=1;
  80     control_value=0.5;
  81     break;
  82   case tBeta:
  83     control_min=-30*degtorad;
  84     control_max=30*degtorad;
  85     control_convert=radtodeg;
  86     break;
  87   case tAlpha:
  88     control_min=fdmex->GetAircraft()->GetAlphaCLMin();
  89     control_max=fdmex->GetAircraft()->GetAlphaCLMax();
  90     if(control_max <= control_min) {
  91       control_max=20*degtorad;
  92       control_min=-5*degtorad;
  93     }
  94     control_value= (control_min+control_max)/2;
  95     control_convert=radtodeg;
  96     solver_eps=tolerance/100;
  97     break;
  98   case tPitchTrim:
  99   case tElevator:
 100   case tRollTrim:
 101   case tAileron:
 102   case tYawTrim:
 103   case tRudder:
 104     control_min=-1;
 105     control_max=1;
 106     state_convert=radtodeg;
 107     solver_eps=tolerance/100;
 108     break;
 109   case tAltAGL:
 110     control_min=0;
 111     control_max=30;
 112     control_value=fdmex->GetPosition()->GetDistanceAGL();
 113     solver_eps=tolerance/100;
 114     break;
 115   case tTheta:
 116     control_min=fdmex->GetRotation()->Gettht() - 5*degtorad;
 117     control_max=fdmex->GetRotation()->Gettht() + 5*degtorad;
 118     state_convert=radtodeg;
 119     break;
 120   case tPhi:
 121     control_min=fdmex->GetRotation()->Getphi() - 30*degtorad;
 122     control_max=fdmex->GetRotation()->Getphi() + 30*degtorad;
 123     state_convert=radtodeg;
 124     control_convert=radtodeg;
 125     break;
 126   case tGamma:
 127     solver_eps=tolerance/100;
 128     control_min=-80*degtorad;
 129     control_max=80*degtorad;
 130     control_convert=radtodeg;
 131     break;
 132   case tHeading:
 133     control_min=fdmex->GetRotation()->Getpsi() - 30*degtorad;
 134     control_max=fdmex->GetRotation()->Getpsi() + 30*degtorad;
 135     state_convert=radtodeg;
 136     break;
 137   }
 138
 139
 140   if (debug_lvl & 2) cout << "Instantiated: FGTrimAxis" << endl;
 141 }
 142
 143 /*****************************************************************************/
 144
 145 FGTrimAxis::~FGTrimAxis()
 146 {
 147   if (debug_lvl & 2) cout << "Destroyed:    FGTrimAxis" << endl;
 148 }
 149
 150 /*****************************************************************************/
 151
 152 void FGTrimAxis::getState(void) {
 153   switch(state) {
 154   case tUdot: state_value=fdmex->GetAircraft()->GetBodyAccel()(1); break;
 155   case tVdot: state_value=fdmex->GetAircraft()->GetBodyAccel()(2); break;
 156   case tWdot: state_value=fdmex->GetAircraft()->GetBodyAccel()(3); break;
 157   case tQdot: state_value=fdmex->GetRotation()->GetPQRdot(2);break;
 158   case tPdot: state_value=fdmex->GetRotation()->GetPQRdot(1); break;
 159   case tRdot: state_value=fdmex->GetRotation()->GetPQRdot(3); break;
 160   case tHmgt: state_value=computeHmgt(); break;
 161   }
 162 }
 163
 164 /*****************************************************************************/
 165
 166 //States are not settable
 167
 168 void FGTrimAxis::getControl(void) {
 169   switch(control) {
 170   case tThrottle:  control_value=fdmex->GetFCS()->GetThrottleCmd(0); break;
 171   case tBeta:      control_value=fdmex->GetTranslation()->Getalpha(); break;
 172   case tAlpha:     control_value=fdmex->GetTranslation()->Getbeta();  break;
 173   case tPitchTrim: control_value=fdmex->GetFCS() -> GetPitchTrimCmd(); break;
 174   case tElevator:  control_value=fdmex->GetFCS() -> GetDeCmd(); break;
 175   case tRollTrim:
 176   case tAileron:   control_value=fdmex->GetFCS() -> GetDaCmd(); break;
 177   case tYawTrim:
 178   case tRudder:    control_value=fdmex->GetFCS() -> GetDrCmd(); break;
 179   case tAltAGL:    control_value=fdmex->GetPosition()->GetDistanceAGL();break;
 180   case tTheta:     control_value=fdmex->GetRotation()->Gettht(); break;
 181   case tPhi:       control_value=fdmex->GetRotation()->Getphi(); break;
 182   case tGamma:     control_value=fdmex->GetPosition()->GetGamma();break;
 183   case tHeading:   control_value=fdmex->GetRotation()->Getpsi(); break;
 184   }
 185 }
 186
 187 /*****************************************************************************/
 188
 189 float FGTrimAxis::computeHmgt(void) {
 190   float diff;
 191
 192   diff   = fdmex->GetRotation()->Getpsi() -
 193              fdmex->GetPosition()->GetGroundTrack();
 194
 195   if( diff < -M_PI ) {
 196      return (diff + 2*M_PI);
 197   } else if( diff > M_PI ) {
 198      return (diff - 2*M_PI);
 199   } else {
 200      return diff;
 201   }
 202
 203 }
 204
 205 /*****************************************************************************/
 206
 207
 208 void FGTrimAxis::setControl(void) {
 209   switch(control) {
 210   case tThrottle:  setThrottlesPct(); break;
 211   case tBeta:      fgic->SetBetaRadIC(control_value); break;
 212   case tAlpha:     fgic->SetAlphaRadIC(control_value);  break;
 213   case tPitchTrim: fdmex->GetFCS()->SetPitchTrimCmd(control_value); break;
 214   case tElevator:  fdmex->GetFCS()->SetDeCmd(control_value); break;
 215   case tRollTrim:
 216   case tAileron:   fdmex->GetFCS()->SetDaCmd(control_value); break;
 217   case tYawTrim:
 218   case tRudder:    fdmex->GetFCS()->SetDrCmd(control_value); break;
 219   case tAltAGL:    fgic->SetAltitudeAGLFtIC(control_value); break;
 220   case tTheta:     fgic->SetPitchAngleRadIC(control_value); break;
 221   case tPhi:       fgic->SetRollAngleRadIC(control_value); break;
 222   case tGamma:     fgic->SetFlightPathAngleRadIC(control_value); break;
 223   case tHeading:   fgic->SetTrueHeadingRadIC(control_value); break;
 224   }
 225 }
 226
 227
 228
 229
 230
 231 /*****************************************************************************/
 232
 233 // the aircraft center of rotation is no longer the cg once the gear
 234 // contact the ground so the altitude needs to be changed when pitch
 235 // and roll angle are adjusted.  Instead of attempting to calculate the
 236 // new center of rotation, pick a gear unit as a reference and use its
 237 // location vector to calculate the new height change. i.e. new altitude =
 238 // earth z component of that vector (which is in body axes )
 239 void FGTrimAxis::SetThetaOnGround(float ff) {
 240   int center,i,ref;
 241
 242   // favor an off-center unit so that the same one can be used for both
 243   // pitch and roll.  An on-center unit is used (for pitch)if that's all
 244   // that's in contact with the ground.
 245   i=0; ref=-1; center=-1;
 246   while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
 247     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) {
 248       if(fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01)
 249         ref=i;
 250       else
 251         center=i;
 252     }
 253     i++;
 254   }
 255   if((ref < 0) && (center >= 0)) {
 256     ref=center;
 257   }
 258   cout << "SetThetaOnGround ref gear: " << ref << endl;
 259   if(ref >= 0) {
 260     float sp=fdmex->GetRotation()->GetSinphi();
 261     float cp=fdmex->GetRotation()->GetCosphi();
 262     float lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
 263     float ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
 264     float lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
 265     float hagl = -1*lx*sin(ff) +
 266                     ly*sp*cos(ff) +
 267                     lz*cp*cos(ff);
 268
 269     fgic->SetAltitudeAGLFtIC(hagl);
 270     cout << "SetThetaOnGround new alt: " << hagl << endl;
 271   }
 272   fgic->SetPitchAngleRadIC(ff);
 273   cout << "SetThetaOnGround new theta: " << ff << endl;
 274 }
 275
 276 /*****************************************************************************/
 277
 278 bool FGTrimAxis::initTheta(void) {
 279   int i,N,iAft, iForward;
 280   float zAft,zForward,zDiff,theta;
 281   bool level;
 282   float saveAlt;
 283
 284   saveAlt=fgic->GetAltitudeAGLFtIC();
 285   fgic->SetAltitudeAGLFtIC(100);
 286
 287
 288   N=fdmex->GetGroundReactions()->GetNumGearUnits();
 289
 290   //find the first wheel unit forward of the cg
 291   //the list is short so a simple linear search is fine
 292   for( i=0; i<N; i++ ) {
 293     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) > 0 ) {
 294         iForward=i;
 295         break;
 296     }
 297   }
 298   //now find the first wheel unit aft of the cg
 299   for( i=0; i<N; i++ ) {
 300     if(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(1) < 0 ) {
 301         iAft=i;
 302         break;
 303     }
 304   }
 305
 306   // now adjust theta till the wheels are the same distance from the ground
 307   zAft=fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear(3);
 308   zForward=fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear(3);
 309   zDiff = zForward - zAft;
 310   level=false;
 311   theta=fgic->GetPitchAngleDegIC();
 312   while(!level && (i < 100)) {
 313         theta+=2.0*zDiff;
 314         fgic->SetPitchAngleDegIC(theta);
 315         fdmex->RunIC(fgic);
 316         zAft=fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear(3);
 317         zForward=fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear(3);
 318         zDiff = zForward - zAft;
 319         //cout << endl << theta << "  " << zDiff << endl;
 320         //cout << "0: " << fdmex->GetGroundReactions()->GetGearUnit(0)->GetLocalGear() << endl;
 321         //cout << "1: " << fdmex->GetGroundReactions()->GetGearUnit(1)->GetLocalGear() << endl;
 322
 323         if(fabs(zDiff ) < 0.1)
 324             level=true;
 325         i++;
 326   }
 327   //cout << i << endl;
 328   cout << "    Initial Theta: " << fdmex->GetRotation()->Gettht()*radtodeg << endl;
 329   control_min=(theta+5)*degtorad;
 330   control_max=(theta-5)*degtorad;
 331   fgic->SetAltitudeAGLFtIC(saveAlt);
 332   if(i < 100)
 333     return true;
 334   else
 335     return false;
 336 }
 337
 338 /*****************************************************************************/
 339
 340 void FGTrimAxis::SetPhiOnGround(float ff) {
 341   int i,ref;
 342
 343   i=0; ref=-1;
 344   //must have an off-center unit here
 345   while( (ref < 0) && (i < fdmex->GetGroundReactions()->GetNumGearUnits()) ) {
 346     if( (fdmex->GetGroundReactions()->GetGearUnit(i)->GetWOW()) &&
 347       (fabs(fdmex->GetGroundReactions()->GetGearUnit(i)->GetBodyLocation(2)) > 0.01))
 348         ref=i;
 349     i++;
 350   }
 351   if(ref >= 0) {
 352     float st=fdmex->GetRotation()->GetSintht();
 353     float ct=fdmex->GetRotation()->GetCostht();
 354     float lx=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(1);
 355     float ly=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(2);
 356     float lz=fdmex->GetGroundReactions()->GetGearUnit(ref)->GetBodyLocation(3);
 357     float hagl = -1*lx*st +
 358                     ly*sin(ff)*ct +
 359                     lz*cos(ff)*ct;
 360
 361     fgic->SetAltitudeAGLFtIC(hagl);
 362   }
 363   fgic->SetRollAngleRadIC(ff);
 364
 365 }
 366
 367 /*****************************************************************************/
 368
 369 void FGTrimAxis::Run(void) {
 370
 371   float last_state_value;
 372   int i;
 373   setControl();
 374   //cout << "FGTrimAxis::Run: " << control_value << endl;
 375   i=0;
 376   bool stable=false;
 377   while(!stable) {
 378     i++;
 379     last_state_value=state_value;
 380     fdmex->RunIC(fgic);
 381     getState();
 382     if(i > 1) {
 383       if((fabs(last_state_value - state_value) < tolerance) || (i >= 100) )
 384         stable=true;
 385     }
 386   }
 387
 388   its_to_stable_value=i;
 389   total_stability_iterations+=its_to_stable_value;
 390   total_iterations++;
 391 }
 392
 393 /*****************************************************************************/
 394
 395 void FGTrimAxis::setThrottlesPct(void) {
 396   float tMin,tMax;
 397   for(unsigned i=0;i<fdmex->GetPropulsion()->GetNumEngines();i++) {
 398       tMin=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMin();
 399       tMax=fdmex->GetPropulsion()->GetEngine(i)->GetThrottleMax();
 400       //cout << "setThrottlespct: " << i << ", " << control_min << ", " << control_max << ", " << control_value;
 401       fdmex->GetFCS()->SetThrottleCmd(i,tMin+control_value*(tMax-tMin));
 402       //cout << "setThrottlespct: " << fdmex->GetFCS()->GetThrottleCmd(i) << endl;
 403       fdmex->RunIC(fgic); //apply throttle change
 404       fdmex->GetPropulsion()->GetSteadyState();
 405   }
 406 }
 407
 408 /*****************************************************************************/
 409
 410 void FGTrimAxis::AxisReport(void) {
 411
 412   char out[80];
 413   sprintf(out,"  %20s: %6.2f %5s: %9.2e Tolerance: %3.0e\n",
 414            GetControlName().c_str(), GetControl()*control_convert,
 415            GetStateName().c_str(), GetState(), GetTolerance());
 416   cout << out;
 417
 418 }
 419
 420 /*****************************************************************************/
 421
 422 float FGTrimAxis::GetAvgStability( void ) {
 423   if(total_iterations > 0) {
 424     return float(total_stability_iterations)/float(total_iterations);
 425   }
 426   return 0;
 427 }
 428
 429 /*****************************************************************************/
 430
 431 void FGTrimAxis::Debug(void)
 432 {
 433 }
 434