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