src/FDM/JSBSim/FGAircraft.cpp

   1 /*******************************************************************************
   2
   3  Module:       FGAircraft.cpp
   4  Author:       Jon S. Berndt
   5  Date started: 12/12/98
   6  Purpose:      Encapsulates an aircraft
   7  Called by:    FGFDMExec
   8
   9  ------------- Copyright (C) 1999  Jon S. Berndt (jsb@hal-pc.org) -------------
  10
  11  This program is free software; you can redistribute it and/or modify it under
  12  the terms of the GNU General Public License as published by the Free Software
  13  Foundation; either version 2 of the License, or (at your option) any later
  14  version.
  15
  16  This program is distributed in the hope that it will be useful, but WITHOUT
  17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  18  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  19  details.
  20
  21  You should have received a copy of the GNU General Public License along with
  22  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  23  Place - Suite 330, Boston, MA  02111-1307, USA.
  24
  25  Further information about the GNU General Public License can also be found on
  26  the world wide web at http://www.gnu.org.
  27
  28 FUNCTIONAL DESCRIPTION
  29 --------------------------------------------------------------------------------
  30 Models the aircraft reactions and forces. This class is instantiated by the
  31 FGFDMExec class and scheduled as an FDM entry. LoadAircraft() is supplied with a
  32 name of a valid, registered aircraft, and the data file is parsed.
  33
  34 HISTORY
  35 --------------------------------------------------------------------------------
  36 12/12/98   JSB   Created
  37 04/03/99   JSB   Changed Aero() method to correct body axis force calculation
  38                  from wind vector. Fix provided by Tony Peden.
  39 05/03/99   JSB   Changed (for the better?) the way configurations are read in.
  40 9/17/99     TP   Combined force and moment functions. Added aero reference
  41                  point to config file. Added calculations for moments due to
  42                                  difference in cg and aero reference point
  43
  44 ********************************************************************************
  45 COMMENTS, REFERENCES,  and NOTES
  46 ********************************************************************************
  47 [1] Cooke, Zyda, Pratt, and McGhee, "NPSNET: Flight Simulation Dynamic Modeling
  48       Using Quaternions", Presence, Vol. 1, No. 4, pp. 404-420  Naval Postgraduate
  49       School, January 1994
  50 [2] D. M. Henderson, "Euler Angles, Quaternions, and Transformation Matrices",
  51       JSC 12960, July 1977
  52 [3] Richard E. McFarland, "A Standard Kinematic Model for Flight Simulation at
  53       NASA-Ames", NASA CR-2497, January 1975
  54 [4] Barnes W. McCormick, "Aerodynamics, Aeronautics, and Flight Mechanics",
  55       Wiley & Sons, 1979 ISBN 0-471-03032-5
  56 [5] Bernard Etkin, "Dynamics of Flight, Stability and Control", Wiley & Sons,
  57       1982 ISBN 0-471-08936-2
  58
  59 The aerodynamic coefficients used in this model are:
  60
  61 Longitudinal
  62   CL0 - Reference lift at zero alpha
  63   CD0 - Reference drag at zero alpha
  64   CDM - Drag due to Mach
  65   CLa - Lift curve slope (w.r.t. alpha)
  66   CDa - Drag curve slope (w.r.t. alpha)
  67   CLq - Lift due to pitch rate
  68   CLM - Lift due to Mach
  69   CLadt - Lift due to alpha rate
  70
  71   Cmadt - Pitching Moment due to alpha rate
  72   Cm0 - Reference Pitching moment at zero alpha
  73   Cma - Pitching moment slope (w.r.t. alpha)
  74   Cmq - Pitch damping (pitch moment due to pitch rate)
  75   CmM - Pitch Moment due to Mach
  76
  77 Lateral
  78   Cyb - Side force due to sideslip
  79   Cyr - Side force due to yaw rate
  80
  81   Clb - Dihedral effect (roll moment due to sideslip)
  82   Clp - Roll damping (roll moment due to roll rate)
  83   Clr - Roll moment due to yaw rate
  84   Cnb - Weathercocking stability (yaw moment due to sideslip)
  85   Cnp - Rudder adverse yaw (yaw moment due to roll rate)
  86   Cnr - Yaw damping (yaw moment due to yaw rate)
  87
  88 Control
  89   CLDe - Lift due to elevator
  90   CDDe - Drag due to elevator
  91   CyDr - Side force due to rudder
  92   CyDa - Side force due to aileron
  93
  94   CmDe - Pitch moment due to elevator
  95   ClDa - Roll moment due to aileron
  96   ClDr - Roll moment due to rudder
  97   CnDr - Yaw moment due to rudder
  98   CnDa - Yaw moment due to aileron
  99
 100 ********************************************************************************
 101 INCLUDES
 102 *******************************************************************************/
 103
 104 #include <sys/stat.h>
 105 #include <sys/types.h>
 106
 107 #ifdef FGFS
 108 # ifndef __BORLANDC__
 109 #  include <Include/compiler.h>
 110 # endif
 111 #  ifdef FG_HAVE_STD_INCLUDES
 112 #    include <cmath>
 113 #  else
 114 #    include <math.h>
 115 #  endif
 116 #else
 117 #  include <cmath>
 118 #endif
 119
 120 #include "FGAircraft.h"
 121 #include "FGTranslation.h"
 122 #include "FGRotation.h"
 123 #include "FGAtmosphere.h"
 124 #include "FGState.h"
 125 #include "FGFDMExec.h"
 126 #include "FGFCS.h"
 127 #include "FGPosition.h"
 128 #include "FGAuxiliary.h"
 129 #include "FGOutput.h"
 130
 131 /*******************************************************************************
 132 ************************************ CODE **************************************
 133 *******************************************************************************/
 134
 135 FGAircraft::FGAircraft(FGFDMExec* fdmex) : FGModel(fdmex)
 136 {
 137   int i;
 138
 139   Name = "FGAircraft";
 140
 141   for (i=0;i<6;i++) coeff_ctr[i] = 0;
 142 }
 143
 144
 145 FGAircraft::~FGAircraft(void)
 146 {
 147 }
 148
 149 bool FGAircraft::LoadAircraft(string aircraft_path, string engine_path, string fname)
 150 {
 151   string path;
 152   string fullpath;
 153   string filename;
 154   string aircraftDef;
 155   string tag;
 156   string holding_string;
 157   char scratch[128];
 158   ifstream coeffInFile;
 159   streampos gpos;
 160   int axis;
 161   string axis_descript;
 162   bool readAeroRp=false;
 163
 164   axis = -1;
 165   aircraftDef = aircraft_path + "/" + fname + "/" + fname + ".cfg";
 166   ifstream aircraftfile(aircraftDef.c_str());
 167   cout << "Reading Aircraft Configuration File: " << aircraftDef << endl;
 168   // Output->SocketStatusOutput("Reading Aircraft Configuration File: " + aircraftDef);
 169
 170   numTanks = numEngines = 0;
 171   numSelectedOxiTanks = numSelectedFuelTanks = 0;
 172
 173   Xcg=Ycg=Zcg=0; //protection for no cg specified in file
 174
 175   while (!aircraftfile.fail()) {
 176     holding_string.erase();
 177     aircraftfile >> holding_string;
 178 #if defined(__BORLANDC__) || defined(FG_HAVE_NATIVE_SGI_COMPILERS) || defined(_MSC_VER)
 179     if (holding_string.compare(0, 2, "//") != 0) {
 180 #else
 181     if (holding_string.compare("//",0,2) != 0) {
 182 #endif
 183       if (holding_string == "CFG_VERSION") {
 184         aircraftfile >> CFGVersion;
 185         cout << "Config file version: " << CFGVersion << endl;
 186         if (CFGVersion < NEEDED_CFG_VERSION) {
 187           cout << endl << "YOU HAVE AN OLD CFG FILE FOR THIS AIRCRAFT."
 188                           " RESULTS WILL BE UNPREDICTABLE !!" << endl << endl;
 189         }
 190       } else if (holding_string == "AIRCRAFT") {
 191         cout << "Reading in Aircraft parameters ..." << endl;
 192       } else if (holding_string == "AERODYNAMICS") {
 193         cout << "Reading in Aerodynamic parameters ..." << endl;
 194       } else if (holding_string == "AC_NAME") {
 195         aircraftfile >> AircraftName;   // String with no embedded spaces
 196         cout << "Aircraft Name: " << AircraftName << endl;
 197       } else if (holding_string == "AC_WINGAREA") {
 198         aircraftfile >> WingArea;
 199         cout << "Aircraft Wing Area: " << WingArea << endl;
 200       } else if (holding_string == "AC_WINGSPAN") {
 201         aircraftfile >> WingSpan;
 202         cout << "Aircraft WingSpan: " << WingSpan << endl;
 203       } else if (holding_string == "AC_CHORD") {
 204         aircraftfile >> cbar;
 205         cout << "Aircraft Chord: " << cbar << endl;
 206       } else if (holding_string == "AC_IXX") {
 207         aircraftfile >> baseIxx;
 208         cout << "Aircraft Base Ixx: " << baseIxx << endl;
 209       } else if (holding_string == "AC_IYY") {
 210         aircraftfile >> baseIyy;
 211         cout << "Aircraft Base Iyy: " << baseIyy << endl;
 212       } else if (holding_string == "AC_IZZ") {
 213         aircraftfile >> baseIzz;
 214         cout << "Aircraft Base Izz: " << baseIzz << endl;
 215       } else if (holding_string == "AC_IXZ") {
 216         aircraftfile >> baseIxz;
 217         cout << "Aircraft Base Ixz: " << baseIxz << endl;
 218       } else if (holding_string == "AC_EMPTYWT") {
 219         aircraftfile >> EmptyWeight;
 220         EmptyMass = EmptyWeight / GRAVITY;
 221         cout << "Aircraft Empty Weight: " << EmptyWeight << endl;
 222       } else if (holding_string == "AC_AERORP") {
 223         aircraftfile >> Xrp >> Yrp >> Zrp;
 224         readAeroRp=true;
 225         cout << "Aerodynamic Reference Point: " << Xrp << " " << Yrp << " " << Zrp << endl;
 226       } else if (holding_string == "AC_CGLOC") {
 227         aircraftfile >> baseXcg >> baseYcg >> baseZcg;
 228         cout << "Aircraft Base C.G.: " << baseXcg << " " << baseYcg << " " << baseZcg << endl;
 229       } else if (holding_string == "AC_EYEPTLOC") {
 230         aircraftfile >> Xep >> Yep >> Zep;
 231         cout << "Pilot Eyepoint: " << Xep << " " << Yep << " " << Zep << endl;
 232       } else if (holding_string == "AC_TANK") {
 233         Tank[numTanks] = new FGTank(aircraftfile);
 234         switch(Tank[numTanks]->GetType()) {
 235         case FGTank::ttFUEL:
 236           numSelectedFuelTanks++;
 237           cout << "Reading in Fuel Tank #" << numSelectedFuelTanks << " parameters ..." << endl;
 238           break;
 239         case FGTank::ttOXIDIZER:
 240           numSelectedOxiTanks++;
 241           cout << "Reading in Oxidizer Tank #" << numSelectedOxiTanks << " parameters ..." << endl;
 242           break;
 243         }
 244         numTanks++;
 245
 246       } else if (holding_string == "AC_GEAR") {
 247
 248         lGear.push_back(new FGLGear(aircraftfile));
 249
 250       } else if (holding_string == "AC_ENGINE") {
 251
 252         aircraftfile >> tag;
 253         cout << "Reading in " << tag << " Engine parameters ..." << endl;
 254         Engine[numEngines] = new FGEngine(FDMExec, engine_path, tag, numEngines);
 255         numEngines++;
 256
 257       } else if (holding_string == "}") {
 258
 259       } else if (holding_string == "{") {
 260
 261       } else if (holding_string == "LIFT") {
 262
 263         axis_descript = "   Lift Coefficients ...";
 264         axis = LiftCoeff;
 265
 266       } else if (holding_string == "DRAG") {
 267
 268         axis_descript = "   Drag Coefficients ...";
 269         axis = DragCoeff;
 270
 271       } else if (holding_string == "SIDE") {
 272
 273         axis_descript = "   Side Coefficients ...";
 274         axis = SideCoeff;
 275
 276       } else if (holding_string == "ROLL") {
 277
 278         axis_descript = "   Roll Coefficients ...";
 279         axis = RollCoeff;
 280
 281       } else if (holding_string == "PITCH") {
 282
 283         axis_descript = "   Pitch Coefficients ...";
 284         axis = PitchCoeff;
 285
 286       } else if (holding_string == "YAW") {
 287
 288         axis_descript = "   Yaw Coefficients ...";
 289         axis = YawCoeff;
 290
 291       }
 292
 293       if (axis >= 0) {
 294         cout << axis_descript << endl;
 295         aircraftfile >> tag;
 296         gpos = aircraftfile.tellg();
 297         aircraftfile >> tag;
 298         if ( !(tag == "}") ) {
 299           while ( !(tag == "}") ) {
 300             aircraftfile.seekg(gpos);
 301             Coeff[axis][coeff_ctr[axis]] = new FGCoefficient(FDMExec, aircraftfile);
 302             coeff_ctr[axis]++;
 303             aircraftfile >> tag;
 304             gpos = aircraftfile.tellg();
 305             aircraftfile >> tag;
 306           }
 307         } else {
 308           cout << "      None found ..." << endl;
 309         }
 310       }
 311       axis = -1;
 312
 313     } else {
 314       aircraftfile.getline(scratch, 127);
 315     }
 316   }
 317
 318   if (!readAeroRp) {
 319     Xrp = Xcg;
 320     Yrp = Ycg;
 321     Zrp = Zcg;
 322     cout << "Aerodynamic reference point not found, set to empty weight cg location" << endl;
 323   }
 324
 325   cout << "End of Configuration File Parsing." << endl;
 326   return true;
 327 }
 328
 329
 330 bool FGAircraft::Run(void)
 331 {
 332   if (!FGModel::Run()) {                 // if false then execute this Run()
 333     GetState();
 334
 335     for (int i = 0; i < 3; i++)  Forces[i] = Moments[i] = 0.0;
 336
 337     MassChange();
 338
 339     FMProp(); FMAero(); FMGear(); FMMass();
 340
 341     PutState();
 342   } else {                               // skip Run() execution this time
 343   }
 344   return false;
 345 }
 346
 347
 348 void FGAircraft::MassChange()
 349 {
 350   float Xt, Xw, Yt, Yw, Zt, Zw, Tw;
 351   float IXXt, IYYt, IZZt, IXZt;
 352   int t;
 353
 354   // UPDATE TANK CONTENTS
 355   //
 356   // For each engine, cycle through the tanks and draw an equal amount of
 357   // fuel (or oxidizer) from each active tank. The needed amount of fuel is
 358   // determined by the engine in the FGEngine class. If more fuel is needed
 359   // than is available in the tank, then that amount is considered a shortage,
 360   // and will be drawn from the next tank. If the engine cannot be fed what it
 361   // needs, it will be considered to be starved, and will shut down.
 362
 363   float Oshortage, Fshortage;
 364
 365   for (int e=0; e<numEngines; e++) {
 366     Fshortage = Oshortage = 0.0;
 367     for (t=0; t<numTanks; t++) {
 368       switch(Engine[e]->GetType()) {
 369       case FGEngine::etRocket:
 370
 371         switch(Tank[t]->GetType()) {
 372         case FGTank::ttFUEL:
 373           if (Tank[t]->GetSelected()) {
 374             Fshortage = Tank[t]->Reduce((Engine[e]->CalcFuelNeed()/
 375                                          numSelectedFuelTanks)*(dt*rate) + Fshortage);
 376           }
 377           break;
 378         case FGTank::ttOXIDIZER:
 379           if (Tank[t]->GetSelected()) {
 380             Oshortage = Tank[t]->Reduce((Engine[e]->CalcOxidizerNeed()/
 381                                          numSelectedOxiTanks)*(dt*rate) + Oshortage);
 382           }
 383           break;
 384         }
 385         break;
 386
 387       case FGEngine::etPiston:
 388       case FGEngine::etTurboJet:
 389       case FGEngine::etTurboProp:
 390
 391         if (Tank[t]->GetSelected()) {
 392           Fshortage = Tank[t]->Reduce((Engine[e]->CalcFuelNeed()/
 393                                        numSelectedFuelTanks)*(dt*rate) + Fshortage);
 394         }
 395         break;
 396       }
 397     }
 398     if ((Fshortage <= 0.0) || (Oshortage <= 0.0)) Engine[e]->SetStarved();
 399     else Engine[e]->SetStarved(false);
 400   }
 401
 402   Weight = EmptyWeight;
 403   for (t=0; t<numTanks; t++)
 404     Weight += Tank[t]->GetContents();
 405
 406   Mass = Weight / GRAVITY;
 407
 408   // Calculate new CG here.
 409
 410   Xt = Yt = Zt = Tw = 0;
 411   Xw = Yw = Zw = 0;
 412   for (t=0; t<numTanks; t++) {
 413     Xt += Tank[t]->GetX()*Tank[t]->GetContents();
 414     Yt += Tank[t]->GetY()*Tank[t]->GetContents();
 415     Zt += Tank[t]->GetZ()*Tank[t]->GetContents();
 416
 417     Tw += Tank[t]->GetContents();
 418   }
 419
 420   Xcg = (Xt + EmptyWeight*baseXcg) / (Tw + EmptyWeight);
 421   Ycg = (Yt + EmptyWeight*baseYcg) / (Tw + EmptyWeight);
 422   Zcg = (Zt + EmptyWeight*baseZcg) / (Tw + EmptyWeight);
 423
 424   // Calculate new moments of inertia here
 425
 426   IXXt = IYYt = IZZt = IXZt = 0.0;
 427   for (t=0; t<numTanks; t++) {
 428     IXXt += ((Tank[t]->GetX()-Xcg)/12.0)*((Tank[t]->GetX() - Xcg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 429     IYYt += ((Tank[t]->GetY()-Ycg)/12.0)*((Tank[t]->GetY() - Ycg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 430     IZZt += ((Tank[t]->GetZ()-Zcg)/12.0)*((Tank[t]->GetZ() - Zcg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 431     IXZt += ((Tank[t]->GetX()-Xcg)/12.0)*((Tank[t]->GetZ() - Zcg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 432   }
 433
 434   Ixx = baseIxx + IXXt;
 435   Iyy = baseIyy + IYYt;
 436   Izz = baseIzz + IZZt;
 437   Ixz = baseIxz + IXZt;
 438
 439 }
 440
 441
 442 void FGAircraft::FMAero(void)
 443 {
 444   float F[3];
 445   float dxcg,dycg,dzcg;
 446   float ca, cb, sa, sb;
 447   int axis_ctr,ctr;
 448   F[0] = F[1] = F[2] = 0.0;
 449
 450   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 451     for (ctr=0; ctr < coeff_ctr[axis_ctr]; ctr++) {
 452       F[axis_ctr] += Coeff[axis_ctr][ctr]->TotalValue();
 453 //      Coeff[axis_ctr][ctr]->DumpSD();
 454     }
 455   }
 456
 457   ca = cos(alpha);
 458   sa = sin(alpha);
 459   cb = cos(beta);
 460   sb = sin(beta);
 461
 462   Forces[0] += - F[DragCoeff]*ca*cb
 463                - F[SideCoeff]*ca*sb
 464                + F[LiftCoeff]*sa;
 465   Forces[1] +=   F[DragCoeff]*sb
 466                  + F[SideCoeff]*cb;
 467   Forces[2] += - F[DragCoeff]*sa*cb
 468                - F[SideCoeff]*sa*sb
 469                - F[LiftCoeff]*ca;
 470
 471   // The d*cg distances below, given in inches, are the distances FROM the c.g.
 472   // TO the reference point. Since the c.g. and ref point are given in inches in
 473   // the structural system (X positive rearwards) and the body coordinate system
 474   // is given with X positive out the nose, the dxcg and dzcg values are
 475   // *rotated* 180 degrees about the Y axis.
 476
 477   dxcg = -(Xrp - Xcg)/12; //cg and rp values are in inches
 478   dycg =  (Yrp - Ycg)/12;
 479   dzcg = -(Zrp - Zcg)/12;
 480
 481   Moments[0] +=  Forces[2]*dycg - Forces[1]*dzcg; //rolling moment
 482   Moments[1] +=  Forces[0]*dzcg - Forces[2]*dxcg; //pitching moment
 483   Moments[2] += -Forces[0]*dycg + Forces[1]*dxcg; //yawing moment
 484
 485   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 486     for (ctr = 0; ctr < coeff_ctr[axis_ctr+3]; ctr++) {
 487       Moments[axis_ctr] += Coeff[axis_ctr+3][ctr]->TotalValue();
 488 //      Coeff[axis_ctr+3][ctr]->DumpSD();
 489     }
 490   }
 491 }
 492
 493
 494 void FGAircraft::FMGear(void)
 495 {
 496   if (GearUp) {
 497     // crash routine
 498   } else {
 499     for (int i=0;i<lGear.size();i++) {
 500       //      lGear[i].
 501     }
 502   }
 503 }
 504
 505
 506 void FGAircraft::FMMass(void)
 507 {
 508   Forces[0] += -GRAVITY*sin(tht) * Mass;
 509   Forces[1] +=  GRAVITY*sin(phi)*cos(tht) * Mass;
 510   Forces[2] +=  GRAVITY*cos(phi)*cos(tht) * Mass;
 511 }
 512
 513
 514 void FGAircraft::FMProp(void)
 515 {
 516   for (int i=0;i<numEngines;i++) {
 517     Forces[0] += Engine[i]->CalcThrust();
 518   }
 519 }
 520
 521 void FGAircraft::GetState(void)
 522 {
 523   dt = State->Getdt();
 524
 525   alpha = Translation->Getalpha();
 526   beta = Translation->Getbeta();
 527   phi = Rotation->Getphi();
 528   tht = Rotation->Gettht();
 529   psi = Rotation->Getpsi();
 530 }
 531
 532
 533 void FGAircraft::PutState(void)
 534 {
 535 }
 536