JSBsim/FGCoefficient.cpp

   1 /*******************************************************************************
   2
   3  Module:       FGCoefficient.cpp
   4  Author:       Jon S. Berndt
   5  Date started: 12/28/98
   6  Purpose:      Encapsulates the stability derivative class FGCoefficient;
   7  Called by:    FGAircraft
   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 This class models the stability derivative coefficient lookup tables or
  31 equations. Note that the coefficients need not be calculated each delta-t.
  32
  33 The coefficient files are located in the axis subdirectory for each aircraft.
  34 For instance, for the X-15, you would find subdirectories under the
  35 aircraft/X-15/ directory named CLIFT, CDRAG, CSIDE, CROLL, CPITCH, CYAW. Under
  36 each of these directories would be files named a, a0, q, and so on. The file
  37 named "a" under the CLIFT directory would contain data for the stability
  38 derivative modeling lift due to a change in alpha. See the FGAircraft.cpp file
  39 for additional information. The coefficient files have the following format:
  40
  41 <name of coefficient>
  42 <short description of coefficient with no embedded spaces>
  43 <method used in calculating the coefficient: TABLE | EQUATION | VECTOR | VALUE>
  44   <parameter identifier for table row (if required)>
  45   <parameter identifier for table column (if required)>
  46 <OR'ed list of parameter identifiers needed to turn this coefficient into a force>
  47 <number of rows in table (if required)>
  48 <number of columns in table (if required)>
  49
  50 <value of parameter indexing into the column of a table or vector - or value
  51   itself for a VALUE coefficient>
  52 <values of parameter indexing into row of a table if TABLE type> <Value of
  53   coefficient at this row and column>
  54
  55 <... repeat above for each column of data in table ...>
  56
  57 As an example for the X-15, for the lift due to mach:
  58
  59 CLa0
  60 Lift_at_zero_alpha
  61 Table 8 3
  62 16384
  63 32768
  64 16387
  65
  66 0.00
  67 0.0 0.0
  68 0.5 0.4
  69 0.9 0.9
  70 1.0 1.6
  71 1.1 1.3
  72 1.4 1.0
  73 2.0 0.5
  74 3.0 0.5
  75
  76 30000.00
  77 0.0 0.0
  78 0.5 0.5
  79 0.9 1.0
  80 1.0 1.7
  81 1.1 1.4
  82 1.4 1.1
  83 2.0 0.6
  84 3.0 0.6
  85
  86 70000.00
  87 0.0 0.0
  88 0.5 0.6
  89 0.9 1.1
  90 1.0 1.7
  91 1.1 1.5
  92 1.4 1.2
  93 2.0 0.7
  94 3.0 0.7
  95
  96 Note that the values in a row which index into the table must be the same value
  97 for each column of data, so the first column of numbers for each altitude are
  98 seen to be equal, and there are the same number of values for each altitude.
  99
 100 See the header file FGCoefficient.h for the values of the identifiers.
 101
 102 HISTORY
 103 --------------------------------------------------------------------------------
 104 12/28/98   JSB   Created
 105
 106 ********************************************************************************
 107 INCLUDES
 108 *******************************************************************************/
 109
 110 #include <stdio.h>
 111 #include <stdlib.h>
 112 #include "FGCoefficient.h"
 113
 114 #include "FGAtmosphere.h"
 115 #include "FGState.h"
 116 #include "FGFDMExec.h"
 117 #include "FGFCS.h"
 118 #include "FGAircraft.h"
 119 #include "FGTranslation.h"
 120 #include "FGRotation.h"
 121 #include "FGPosition.h"
 122 #include "FGAuxiliary.h"
 123 #include "FGOutput.h"
 124
 125 /*******************************************************************************
 126 ************************************ CODE **************************************
 127 *******************************************************************************/
 128
 129 FGCoefficient::FGCoefficient(FGFDMExec* fdex)
 130 {
 131   FDMExec     = fdex;
 132   State       = FDMExec->GetState();
 133   Atmosphere  = FDMExec->GetAtmosphere();
 134   FCS         = FDMExec->GetFCS();
 135   Aircraft    = FDMExec->GetAircraft();
 136   Translation = FDMExec->GetTranslation();
 137   Rotation    = FDMExec->GetRotation();
 138   Position    = FDMExec->GetPosition();
 139   Auxiliary   = FDMExec->GetAuxiliary();
 140   Output      = FDMExec->GetOutput();
 141
 142   rows = columns = 0;
 143 }
 144
 145
 146 FGCoefficient::FGCoefficient(FGFDMExec* fdex, char* fname)
 147 {
 148   int r, c;
 149   float ftrashcan;
 150
 151   FDMExec     = fdex;
 152   State       = FDMExec->GetState();
 153   Atmosphere  = FDMExec->GetAtmosphere();
 154   FCS         = FDMExec->GetFCS();
 155   Aircraft    = FDMExec->GetAircraft();
 156   Translation = FDMExec->GetTranslation();
 157   Rotation    = FDMExec->GetRotation();
 158   Position    = FDMExec->GetPosition();
 159   Auxiliary   = FDMExec->GetAuxiliary();
 160   Output      = FDMExec->GetOutput();
 161
 162   ifstream coeffDefFile(fname);
 163
 164   if (coeffDefFile) {
 165     if (!coeffDefFile.fail()) {
 166       coeffDefFile >> name;
 167       coeffDefFile >> description;
 168       coeffDefFile >> method;
 169
 170       if (strcmp(method,"EQUATION") == 0) type = 4;
 171       else if (strcmp(method,"TABLE") == 0) type = 3;
 172       else if (strcmp(method,"VECTOR") == 0) type = 2;
 173       else if (strcmp(method,"VALUE") == 0) type = 1;
 174       else type = 0;
 175
 176       if (type == 2 || type == 3) {
 177         coeffDefFile >> rows;
 178         if (type == 3) {
 179           coeffDefFile >> columns;
 180         }
 181         coeffDefFile >> LookupR;
 182       }
 183
 184       if (type == 3) {
 185         coeffDefFile >> LookupC;
 186       }
 187
 188       coeffDefFile >> multipliers;
 189
 190       mult_count = 0;
 191       if (multipliers & FG_QBAR) {
 192         mult_idx[mult_count] = FG_QBAR;
 193         mult_count++;
 194       }
 195       if (multipliers & FG_WINGAREA) {
 196         mult_idx[mult_count] = FG_WINGAREA;
 197         mult_count++;
 198       }
 199       if (multipliers & FG_WINGSPAN) {
 200         mult_idx[mult_count] = FG_WINGSPAN;
 201         mult_count++;
 202       }
 203       if (multipliers & FG_CBAR) {
 204         mult_idx[mult_count] = FG_CBAR;
 205         mult_count++;
 206       }
 207       if (multipliers & FG_ALPHA) {
 208         mult_idx[mult_count] = FG_ALPHA;
 209         mult_count++;
 210       }
 211       if (multipliers & FG_ALPHADOT) {
 212         mult_idx[mult_count] = FG_ALPHADOT;
 213         mult_count++;
 214       }
 215       if (multipliers & FG_BETA) {
 216         mult_idx[mult_count] = FG_BETA;
 217         mult_count++;
 218       }
 219       if (multipliers & FG_BETADOT) {
 220         mult_idx[mult_count] = FG_BETADOT;
 221         mult_count++;
 222       }
 223       if (multipliers & FG_PITCHRATE) {
 224         mult_idx[mult_count] = FG_PITCHRATE;
 225         mult_count++;
 226       }
 227       if (multipliers & FG_ROLLRATE) {
 228         mult_idx[mult_count] = FG_ROLLRATE;
 229         mult_count++;
 230       }
 231       if (multipliers & FG_YAWRATE) {
 232         mult_idx[mult_count] = FG_YAWRATE;
 233         mult_count++;
 234       }
 235       if (multipliers & FG_ELEVATOR) {
 236         mult_idx[mult_count] = FG_ELEVATOR;
 237         mult_count++;
 238       }
 239       if (multipliers & FG_AILERON) {
 240         mult_idx[mult_count] = FG_AILERON;
 241         mult_count++;
 242       }
 243       if (multipliers & FG_RUDDER) {
 244         mult_idx[mult_count] = FG_RUDDER;
 245         mult_count++;
 246       }
 247       if (multipliers & FG_MACH) {
 248         mult_idx[mult_count] = FG_MACH;
 249         mult_count++;
 250       }
 251       if (multipliers & FG_ALTITUDE) {
 252         mult_idx[mult_count] = FG_ALTITUDE;
 253         mult_count++;
 254       }
 255
 256       switch(type) {
 257       case 1:
 258         coeffDefFile >> StaticValue;
 259         break;
 260       case 2:
 261         Allocate(rows,2);
 262
 263         for (r=1;r<=rows;r++) {
 264           coeffDefFile >> Table3D[r][0];
 265           coeffDefFile >> Table3D[r][1];
 266         }
 267         break;
 268       case 3:
 269         Allocate(rows, columns);
 270
 271         for (c=1;c<=columns;c++) {
 272           coeffDefFile >> Table3D[0][c];
 273           for (r=1;r<=rows;r++) {
 274             if ( c==1 ) coeffDefFile >> Table3D[r][0];
 275             else coeffDefFile >> ftrashcan;
 276             coeffDefFile >> Table3D[r][c];
 277           }
 278         }
 279         break;
 280       }
 281     } else {
 282       cerr << "Empty file" << endl;
 283     }
 284     coeffDefFile.close();
 285   }
 286 }
 287
 288
 289 FGCoefficient::FGCoefficient(FGFDMExec* fdex, int r, int c)
 290 {
 291   FDMExec     = fdex;
 292   State       = FDMExec->GetState();
 293   Atmosphere  = FDMExec->GetAtmosphere();
 294   FCS         = FDMExec->GetFCS();
 295   Aircraft    = FDMExec->GetAircraft();
 296   Translation = FDMExec->GetTranslation();
 297   Rotation    = FDMExec->GetRotation();
 298   Position    = FDMExec->GetPosition();
 299   Auxiliary   = FDMExec->GetAuxiliary();
 300   Output      = FDMExec->GetOutput();
 301
 302   rows = r;
 303   columns = c;
 304   Allocate(r,c);
 305 }
 306
 307
 308 FGCoefficient::FGCoefficient(FGFDMExec* fdex, int n)
 309 {
 310   FDMExec     = fdex;
 311   State       = FDMExec->GetState();
 312   Atmosphere  = FDMExec->GetAtmosphere();
 313   FCS         = FDMExec->GetFCS();
 314   Aircraft    = FDMExec->GetAircraft();
 315   Translation = FDMExec->GetTranslation();
 316   Rotation    = FDMExec->GetRotation();
 317   Position    = FDMExec->GetPosition();
 318   Auxiliary   = FDMExec->GetAuxiliary();
 319   Output      = FDMExec->GetOutput();
 320
 321   rows = n;
 322   columns = 0;
 323   Allocate(n);
 324 }
 325
 326
 327 FGCoefficient::~FGCoefficient(void)
 328 {
 329 }
 330
 331
 332 bool FGCoefficient::Allocate(int r, int c)
 333 {
 334   rows = r;
 335   columns = c;
 336   Table3D = new float*[r+1];
 337   for (int i=0;i<=r;i++) Table3D[i] = new float[c+1];
 338   return true;
 339 }
 340
 341
 342 bool FGCoefficient::Allocate(int n)
 343 {
 344   rows = n;
 345   columns = 0;
 346   Table2D = new float[n+1];
 347   return true;
 348 }
 349
 350
 351 float FGCoefficient::Value(float rVal, float cVal)
 352 {
 353   float rFactor, cFactor, col1temp, col2temp, Value;
 354   int r, c, midx;
 355
 356   if (rows < 2 || columns < 2) return 0.0;
 357
 358   for (r=1;r<=rows;r++) if (Table3D[r][0] >= rVal) break;
 359   for (c=1;c<=columns;c++) if (Table3D[0][c] >= cVal) break;
 360
 361   c = c < 2 ? 2 : (c > columns ? columns : c);
 362   r = r < 2 ? 2 : (r > rows    ? rows    : r);
 363
 364   rFactor = (rVal - Table3D[r-1][0]) / (Table3D[r][0] - Table3D[r-1][0]);
 365   cFactor = (cVal - Table3D[0][c-1]) / (Table3D[0][c] - Table3D[0][c-1]);
 366
 367   col1temp = rFactor*(Table3D[r][c-1] - Table3D[r-1][c-1]) + Table3D[r-1][c-1];
 368   col2temp = rFactor*(Table3D[r][c] - Table3D[r-1][c]) + Table3D[r-1][c];
 369
 370   Value = col1temp + cFactor*(col2temp - col1temp);
 371
 372   for (midx=0;midx<mult_count;midx++) {
 373     Value *= GetCoeffVal(mult_idx[midx]);
 374   }
 375
 376   return Value;
 377 }
 378
 379
 380 float FGCoefficient::Value(float Val)
 381 {
 382   float Factor, Value;
 383   int r, midx;
 384
 385   if (rows < 2) return 0.0;
 386
 387   for (r=1;r<=rows;r++) if (Table3D[r][0] >= Val) break;
 388   r = r < 2 ? 2 : (r > rows    ? rows    : r);
 389
 390   // make sure denominator below does not go to zero.
 391   if (Table3D[r][0] != Table3D[r-1][0]) {
 392     Factor = (Val - Table3D[r-1][0]) / (Table3D[r][0] - Table3D[r-1][0]);
 393   } else {
 394     Factor = 1.0;
 395   }
 396   Value = Factor*(Table3D[r][1] - Table3D[r-1][1]) + Table3D[r-1][1];
 397
 398   for (midx=0;midx<mult_count;midx++) {
 399     Value *= GetCoeffVal(mult_idx[midx]);
 400   }
 401
 402   return Value;
 403 }
 404
 405
 406 float FGCoefficient::Value()
 407 {
 408   switch(type) {
 409   case 0:
 410     return -1;
 411   case 1:
 412     return (StaticValue);
 413   case 2:
 414     return (Value(GetCoeffVal(LookupR)));
 415   case 3:
 416     return (Value(GetCoeffVal(LookupR),GetCoeffVal(LookupC)));
 417   case 4:
 418     return 0.0;
 419   }
 420   return 0;
 421 }
 422
 423 float FGCoefficient::GetCoeffVal(int val_idx)
 424 {
 425   switch(val_idx) {
 426   case FG_QBAR:
 427     return State->Getqbar();
 428   case FG_WINGAREA:
 429     return Aircraft->GetWingArea();
 430   case FG_WINGSPAN:
 431     return Aircraft->GetWingSpan();
 432   case FG_CBAR:
 433     return Aircraft->Getcbar();
 434   case FG_ALPHA:
 435     return Translation->Getalpha();
 436   case FG_ALPHADOT:
 437     return State->Getadot();
 438   case FG_BETA:
 439     return Translation->Getbeta();
 440   case FG_BETADOT:
 441     return State->Getbdot();
 442   case FG_PITCHRATE:
 443     return Rotation->GetQ();
 444   case FG_ROLLRATE:
 445     return Rotation->GetP();
 446   case FG_YAWRATE:
 447     return Rotation->GetR();
 448   case FG_ELEVATOR:
 449     return FCS->GetDe();
 450   case FG_AILERON:
 451     return FCS->GetDa();
 452   case FG_RUDDER:
 453     return FCS->GetDr();
 454   case FG_MACH:
 455     return State->GetMach();
 456   case FG_ALTITUDE:
 457     return State->Geth();
 458   }
 459   return 0;
 460 }