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