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/types.h>
 105 #include <sys/stat.h>
 106
 107 #ifdef FGFS
 108 # ifndef __BORLANDC__
 109 #  include <simgear/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(0);
 160   int axis;
 161   string axis_descript;
 162   bool readAeroRp=false;
 163
 164   axis = -1;
 165 #ifdef MACOS
 166   aircraftDef = aircraft_path + ":" + fname + ":" + fname + ".cfg";
 167 #else
 168   aircraftDef = aircraft_path + "/" + fname + "/" + fname + ".cfg";
 169 #endif
 170   ifstream aircraftfile(aircraftDef.c_str());
 171   cout << "Reading Aircraft Configuration File: " << aircraftDef << endl;
 172   // Output->SocketStatusOutput("Reading Aircraft Configuration File: " + aircraftDef);
 173
 174   numTanks = numEngines = 0;
 175   numSelectedOxiTanks = numSelectedFuelTanks = 0;
 176
 177   Xcg=Ycg=Zcg=0; //protection for no cg specified in file
 178
 179   while (!aircraftfile.fail()) {
 180     holding_string.erase();
 181     aircraftfile >> holding_string;
 182 #if defined(__BORLANDC__) || defined(FG_HAVE_NATIVE_SGI_COMPILERS) || defined(_MSC_VER) || defined(__MWERKS__)
 183     if (holding_string.compare(0, 2, "//") != 0) {
 184 #else
 185     if (holding_string.compare("//",0,2) != 0) {
 186 #endif
 187       if (holding_string == "CFG_VERSION") {
 188         aircraftfile >> CFGVersion;
 189         cout << "Config file version: " << CFGVersion << endl;
 190         if (CFGVersion < NEEDED_CFG_VERSION) {
 191           cout << endl << "YOU HAVE AN OLD CFG FILE FOR THIS AIRCRAFT."
 192                           " RESULTS WILL BE UNPREDICTABLE !!" << endl << endl;
 193         }
 194       } else if (holding_string == "AIRCRAFT") {
 195         cout << "Reading in Aircraft parameters ..." << endl;
 196       } else if (holding_string == "AERODYNAMICS") {
 197         cout << "Reading in Aerodynamic parameters ..." << endl;
 198       } else if (holding_string == "AC_NAME") {
 199         aircraftfile >> AircraftName;   // String with no embedded spaces
 200         cout << "Aircraft Name: " << AircraftName << endl;
 201       } else if (holding_string == "AC_WINGAREA") {
 202         aircraftfile >> WingArea;
 203         cout << "Aircraft Wing Area: " << WingArea << endl;
 204       } else if (holding_string == "AC_WINGSPAN") {
 205         aircraftfile >> WingSpan;
 206         cout << "Aircraft WingSpan: " << WingSpan << endl;
 207       } else if (holding_string == "AC_CHORD") {
 208         aircraftfile >> cbar;
 209         cout << "Aircraft Chord: " << cbar << endl;
 210       } else if (holding_string == "AC_IXX") {
 211         aircraftfile >> baseIxx;
 212         cout << "Aircraft Base Ixx: " << baseIxx << endl;
 213       } else if (holding_string == "AC_IYY") {
 214         aircraftfile >> baseIyy;
 215         cout << "Aircraft Base Iyy: " << baseIyy << endl;
 216       } else if (holding_string == "AC_IZZ") {
 217         aircraftfile >> baseIzz;
 218         cout << "Aircraft Base Izz: " << baseIzz << endl;
 219       } else if (holding_string == "AC_IXZ") {
 220         aircraftfile >> baseIxz;
 221         cout << "Aircraft Base Ixz: " << baseIxz << endl;
 222       } else if (holding_string == "AC_EMPTYWT") {
 223         aircraftfile >> EmptyWeight;
 224         EmptyMass = EmptyWeight / GRAVITY;
 225         cout << "Aircraft Empty Weight: " << EmptyWeight << endl;
 226       } else if (holding_string == "AC_AERORP") {
 227         aircraftfile >> Xrp >> Yrp >> Zrp;
 228         readAeroRp=true;
 229         cout << "Aerodynamic Reference Point: " << Xrp << " " << Yrp << " " << Zrp << endl;
 230       } else if (holding_string == "AC_CGLOC") {
 231         aircraftfile >> baseXcg >> baseYcg >> baseZcg;
 232         cout << "Aircraft Base C.G.: " << baseXcg << " " << baseYcg << " " << baseZcg << endl;
 233       } else if (holding_string == "AC_EYEPTLOC") {
 234         aircraftfile >> Xep >> Yep >> Zep;
 235         cout << "Pilot Eyepoint: " << Xep << " " << Yep << " " << Zep << endl;
 236       } else if (holding_string == "AC_TANK") {
 237         Tank[numTanks] = new FGTank(aircraftfile);
 238         switch(Tank[numTanks]->GetType()) {
 239         case FGTank::ttFUEL:
 240           numSelectedFuelTanks++;
 241           cout << "Reading in Fuel Tank #" << numSelectedFuelTanks << " parameters ..." << endl;
 242           break;
 243         case FGTank::ttOXIDIZER:
 244           numSelectedOxiTanks++;
 245           cout << "Reading in Oxidizer Tank #" << numSelectedOxiTanks << " parameters ..." << endl;
 246           break;
 247         }
 248         numTanks++;
 249
 250       } else if (holding_string == "AC_GEAR") {
 251
 252         lGear.push_back(new FGLGear(aircraftfile));
 253
 254       } else if (holding_string == "AC_ENGINE") {
 255
 256         aircraftfile >> tag;
 257         cout << "Reading in " << tag << " Engine parameters ..." << endl;
 258         Engine[numEngines] = new FGEngine(FDMExec, engine_path, tag, numEngines);
 259         numEngines++;
 260
 261       } else if (holding_string == "}") {
 262
 263       } else if (holding_string == "{") {
 264
 265       } else if (holding_string == "LIFT") {
 266
 267         axis_descript = "   Lift Coefficients ...";
 268         axis = LiftCoeff;
 269
 270       } else if (holding_string == "DRAG") {
 271
 272         axis_descript = "   Drag Coefficients ...";
 273         axis = DragCoeff;
 274
 275       } else if (holding_string == "SIDE") {
 276
 277         axis_descript = "   Side Coefficients ...";
 278         axis = SideCoeff;
 279
 280       } else if (holding_string == "ROLL") {
 281
 282         axis_descript = "   Roll Coefficients ...";
 283         axis = RollCoeff;
 284
 285       } else if (holding_string == "PITCH") {
 286
 287         axis_descript = "   Pitch Coefficients ...";
 288         axis = PitchCoeff;
 289
 290       } else if (holding_string == "YAW") {
 291
 292         axis_descript = "   Yaw Coefficients ...";
 293         axis = YawCoeff;
 294
 295       }
 296
 297       if (axis >= 0) {
 298         cout << axis_descript << endl;
 299         aircraftfile >> tag;
 300         gpos = aircraftfile.tellg();
 301         aircraftfile >> tag;
 302         if ( !(tag == "}") ) {
 303           while ( !(tag == "}") ) {
 304             aircraftfile.seekg(gpos);
 305             Coeff[axis][coeff_ctr[axis]] = new FGCoefficient(FDMExec, aircraftfile);
 306             coeff_ctr[axis]++;
 307             aircraftfile >> tag;
 308             gpos = aircraftfile.tellg();
 309             aircraftfile >> tag;
 310           }
 311         } else {
 312           cout << "      None found ..." << endl;
 313         }
 314       }
 315       axis = -1;
 316
 317     } else {
 318       aircraftfile.getline(scratch, 127);
 319     }
 320   }
 321
 322   if (!readAeroRp) {
 323     Xrp = Xcg;
 324     Yrp = Ycg;
 325     Zrp = Zcg;
 326     cout << "Aerodynamic reference point not found, set to empty weight cg location" << endl;
 327   }
 328
 329   cout << "End of Configuration File Parsing." << endl;
 330   return true;
 331 }
 332
 333
 334 bool FGAircraft::Run(void)
 335 {
 336   if (!FGModel::Run()) {                 // if false then execute this Run()
 337     GetState();
 338
 339     for (int i = 0; i < 3; i++)  Forces[i] = Moments[i] = 0.0;
 340
 341     MassChange();
 342
 343     FMProp(); FMAero(); FMGear(); FMMass();
 344
 345     PutState();
 346   } else {                               // skip Run() execution this time
 347   }
 348   return false;
 349 }
 350
 351
 352 void FGAircraft::MassChange()
 353 {
 354   float Xt, Xw, Yt, Yw, Zt, Zw, Tw;
 355   float IXXt, IYYt, IZZt, IXZt;
 356   int t;
 357
 358   // UPDATE TANK CONTENTS
 359   //
 360   // For each engine, cycle through the tanks and draw an equal amount of
 361   // fuel (or oxidizer) from each active tank. The needed amount of fuel is
 362   // determined by the engine in the FGEngine class. If more fuel is needed
 363   // than is available in the tank, then that amount is considered a shortage,
 364   // and will be drawn from the next tank. If the engine cannot be fed what it
 365   // needs, it will be considered to be starved, and will shut down.
 366
 367   float Oshortage, Fshortage;
 368
 369   for (int e=0; e<numEngines; e++) {
 370     Fshortage = Oshortage = 0.0;
 371     for (t=0; t<numTanks; t++) {
 372       switch(Engine[e]->GetType()) {
 373       case FGEngine::etRocket:
 374
 375         switch(Tank[t]->GetType()) {
 376         case FGTank::ttFUEL:
 377           if (Tank[t]->GetSelected()) {
 378             Fshortage = Tank[t]->Reduce((Engine[e]->CalcFuelNeed()/
 379                                          numSelectedFuelTanks)*(dt*rate) + Fshortage);
 380           }
 381           break;
 382         case FGTank::ttOXIDIZER:
 383           if (Tank[t]->GetSelected()) {
 384             Oshortage = Tank[t]->Reduce((Engine[e]->CalcOxidizerNeed()/
 385                                          numSelectedOxiTanks)*(dt*rate) + Oshortage);
 386           }
 387           break;
 388         }
 389         break;
 390
 391       case FGEngine::etPiston:
 392       case FGEngine::etTurboJet:
 393       case FGEngine::etTurboProp:
 394
 395         if (Tank[t]->GetSelected()) {
 396           Fshortage = Tank[t]->Reduce((Engine[e]->CalcFuelNeed()/
 397                                        numSelectedFuelTanks)*(dt*rate) + Fshortage);
 398         }
 399         break;
 400       }
 401     }
 402     if ((Fshortage <= 0.0) || (Oshortage <= 0.0)) Engine[e]->SetStarved();
 403     else Engine[e]->SetStarved(false);
 404   }
 405
 406   Weight = EmptyWeight;
 407   for (t=0; t<numTanks; t++)
 408     Weight += Tank[t]->GetContents();
 409
 410   Mass = Weight / GRAVITY;
 411
 412   // Calculate new CG here.
 413
 414   Xt = Yt = Zt = Tw = 0;
 415   Xw = Yw = Zw = 0;
 416   for (t=0; t<numTanks; t++) {
 417     Xt += Tank[t]->GetX()*Tank[t]->GetContents();
 418     Yt += Tank[t]->GetY()*Tank[t]->GetContents();
 419     Zt += Tank[t]->GetZ()*Tank[t]->GetContents();
 420
 421     Tw += Tank[t]->GetContents();
 422   }
 423
 424   Xcg = (Xt + EmptyWeight*baseXcg) / (Tw + EmptyWeight);
 425   Ycg = (Yt + EmptyWeight*baseYcg) / (Tw + EmptyWeight);
 426   Zcg = (Zt + EmptyWeight*baseZcg) / (Tw + EmptyWeight);
 427
 428   // Calculate new moments of inertia here
 429
 430   IXXt = IYYt = IZZt = IXZt = 0.0;
 431   for (t=0; t<numTanks; t++) {
 432     IXXt += ((Tank[t]->GetX()-Xcg)/12.0)*((Tank[t]->GetX() - Xcg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 433     IYYt += ((Tank[t]->GetY()-Ycg)/12.0)*((Tank[t]->GetY() - Ycg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 434     IZZt += ((Tank[t]->GetZ()-Zcg)/12.0)*((Tank[t]->GetZ() - Zcg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 435     IXZt += ((Tank[t]->GetX()-Xcg)/12.0)*((Tank[t]->GetZ() - Zcg)/12.0)*Tank[t]->GetContents()/GRAVITY;
 436   }
 437
 438   Ixx = baseIxx + IXXt;
 439   Iyy = baseIyy + IYYt;
 440   Izz = baseIzz + IZZt;
 441   Ixz = baseIxz + IXZt;
 442
 443 }
 444
 445
 446 void FGAircraft::FMAero(void)
 447 {
 448   float F[3];
 449   float dxcg,dycg,dzcg;
 450   float ca, cb, sa, sb;
 451   int axis_ctr,ctr;
 452   F[0] = F[1] = F[2] = 0.0;
 453
 454   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 455     for (ctr=0; ctr < coeff_ctr[axis_ctr]; ctr++) {
 456       F[axis_ctr] += Coeff[axis_ctr][ctr]->TotalValue();
 457 //      Coeff[axis_ctr][ctr]->DumpSD();
 458     }
 459   }
 460
 461   ca = cos(alpha);
 462   sa = sin(alpha);
 463   cb = cos(beta);
 464   sb = sin(beta);
 465
 466   Forces[0] += - F[DragCoeff]*ca*cb
 467                - F[SideCoeff]*ca*sb
 468                + F[LiftCoeff]*sa;
 469   Forces[1] +=   F[DragCoeff]*sb
 470                  + F[SideCoeff]*cb;
 471   Forces[2] += - F[DragCoeff]*sa*cb
 472                - F[SideCoeff]*sa*sb
 473                - F[LiftCoeff]*ca;
 474
 475   // The d*cg distances below, given in inches, are the distances FROM the c.g.
 476   // TO the reference point. Since the c.g. and ref point are given in inches in
 477   // the structural system (X positive rearwards) and the body coordinate system
 478   // is given with X positive out the nose, the dxcg and dzcg values are
 479   // *rotated* 180 degrees about the Y axis.
 480
 481   dxcg = -(Xrp - Xcg)/12; //cg and rp values are in inches
 482   dycg =  (Yrp - Ycg)/12;
 483   dzcg = -(Zrp - Zcg)/12;
 484
 485   Moments[0] +=  Forces[2]*dycg - Forces[1]*dzcg; //rolling moment
 486   Moments[1] +=  Forces[0]*dzcg - Forces[2]*dxcg; //pitching moment
 487   Moments[2] += -Forces[0]*dycg + Forces[1]*dxcg; //yawing moment
 488
 489   for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
 490     for (ctr = 0; ctr < coeff_ctr[axis_ctr+3]; ctr++) {
 491       Moments[axis_ctr] += Coeff[axis_ctr+3][ctr]->TotalValue();
 492 //      Coeff[axis_ctr+3][ctr]->DumpSD();
 493     }
 494   }
 495 }
 496
 497
 498 void FGAircraft::FMGear(void)
 499 {
 500   if (GearUp) {
 501     // crash routine
 502   } else {
 503     for (int i=0;i<lGear.size();i++) {
 504       //      lGear[i].
 505     }
 506   }
 507 }
 508
 509
 510 void FGAircraft::FMMass(void)
 511 {
 512   Forces[0] += -GRAVITY*sin(tht) * Mass;
 513   Forces[1] +=  GRAVITY*sin(phi)*cos(tht) * Mass;
 514   Forces[2] +=  GRAVITY*cos(phi)*cos(tht) * Mass;
 515 }
 516
 517
 518 void FGAircraft::FMProp(void)
 519 {
 520   for (int i=0;i<numEngines;i++) {
 521     Forces[0] += Engine[i]->CalcThrust();
 522   }
 523 }
 524
 525 void FGAircraft::GetState(void)
 526 {
 527   dt = State->Getdt();
 528
 529   alpha = Translation->Getalpha();
 530   beta = Translation->Getbeta();
 531   phi = Rotation->Getphi();
 532   tht = Rotation->Gettht();
 533   psi = Rotation->Getpsi();
 534 }
 535
 536
 537 void FGAircraft::PutState(void)
 538 {
 539 }
 540