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