fly/flightgear
1
0
Fork 0
Code Activity

Merge /u/jsb1685/flightgear/ branch yasim into next

Browse source 
https://sourceforge.net/p/flightgear/flightgear/merge-requests/89/
This commit is contained in:
James Turner 2017-05-02 15:09:13 +00:00
commit 2747a33e6b
32 changed files with 546 additions and 758 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

225
src/FDM/YASim/Airplane.cpp
View file

@ -33,30 +33,7 @@ const float SOLVE_TWEAK = 0.3226;
Airplane::Airplane()
{
_emptyWeight = 0;
_pilotPos[0] = _pilotPos[1] = _pilotPos[2] = 0;
_wing = 0;
_tail = 0;
_ballast = 0;
_cruiseP = 0;
_cruiseT = 0;
_cruiseSpeed = 0;
_cruiseWeight = 0;
_cruiseGlideAngle = 0;
_approachP = 0;
_approachT = 0;
_approachSpeed = 0;
_approachAoA = 0;
_approachWeight = 0;
_approachGlideAngle = 0;
_dragFactor = 1;
_liftRatio = 1;
_cruiseAoA = 0;
_tailIncidence = 0;
_failureMsg = 0;
_wingsN = 0;
_approachConfig.isApproach = true;
}
Airplane::~Airplane()
@ -82,10 +59,10 @@ Airplane::~Airplane()
}
for(i=0; i<_solveWeights.size(); i++)
delete (SolveWeight*)_solveWeights.get(i);
for(i=0; i<_cruiseControls.size(); i++)
delete (Control*)_cruiseControls.get(i);
for(i=0; i<_approachControls.size(); i++) {
Control* c = (Control*)_approachControls.get(i);
for(i=0; i<_cruiseConfig.controls.size(); i++)
delete (Control*)_cruiseConfig.controls.get(i);
for(i=0; i<_approachConfig.controls.size(); i++) {
Control* c = (Control*)_approachConfig.controls.get(i);
if(c != &_approachElevator)
delete c;
}
@ -148,30 +125,29 @@ void Airplane::updateGearState()
void Airplane::setApproach(float speed, float altitude, float aoa, float fuel, float gla)
{
_approachSpeed = speed;
_approachP = Atmosphere::getStdPressure(altitude);
_approachT = Atmosphere::getStdTemperature(altitude);
_approachAoA = aoa;
_approachFuel = fuel;
_approachGlideAngle = gla;
_approachConfig.speed = speed;
_approachConfig.altitude = altitude;
_approachConfig.state.setupOrientationFromAoa(aoa); // see runConfig()
_approachConfig.aoa = aoa; // not strictly needed, see runConfig()
_approachConfig.fuel = fuel;
_approachConfig.glideAngle = gla;
}
void Airplane::setCruise(float speed, float altitude, float fuel, float gla)
{
_cruiseSpeed = speed;
_cruiseP = Atmosphere::getStdPressure(altitude);
_cruiseT = Atmosphere::getStdTemperature(altitude);
_cruiseAoA = 0;
_cruiseConfig.speed = speed;
_cruiseConfig.altitude = altitude;
_cruiseConfig.aoa = 0;
_tailIncidence = 0;
_cruiseFuel = fuel;
_cruiseGlideAngle = gla;
_cruiseConfig.fuel = fuel;
_cruiseConfig.glideAngle = gla;
}
void Airplane::setElevatorControl(int control)
{
_approachElevator.control = control;
_approachElevator.val = 0;
_approachControls.add(&_approachElevator);
_approachConfig.controls.add(&_approachElevator);
}
void Airplane::addApproachControl(int control, float val)
@ -179,7 +155,7 @@ void Airplane::addApproachControl(int control, float val)
Control* c = new Control();
c->control = control;
c->val = val;
_approachControls.add(c);
_approachConfig.controls.add(c);
}
void Airplane::addCruiseControl(int control, float val)
@ -187,7 +163,7 @@ void Airplane::addCruiseControl(int control, float val)
Control* c = new Control();
c->control = control;
c->val = val;
_cruiseControls.add(c);
_cruiseConfig.controls.add(c);
}
void Airplane::addSolutionWeight(bool approach, int idx, float wgt)
@ -300,31 +276,11 @@ void Airplane::setFuelFraction(float frac)
}
}
void Airplane::setupState(const float aoa, const float speed, const float gla, State* s)
{
float cosAoA = Math::cos(aoa);
float sinAoA = Math::sin(aoa);
s->orient[0] = cosAoA; s->orient[1] = 0; s->orient[2] = sinAoA;
s->orient[3] = 0; s->orient[4] = 1; s->orient[5] = 0;
s->orient[6] = -sinAoA; s->orient[7] = 0; s->orient[8] = cosAoA;
//? what is gla? v[1]=y, v[2]=z? should sin go to v2 instead v1?
s->v[0] = speed*Math::cos(gla); s->v[1] = -speed*Math::sin(gla); s->v[2] = 0;
for(int i=0; i<3; i++)
s->pos[i] = s->rot[i] = s->acc[i] = s->racc[i] = 0;
// Put us 1m above the origin, or else the gravity computation in
// Model goes nuts
s->pos[2] = 1;
}
/**
* @brief add contact point for crash detection
* used to add wingtips and fuselage nose and tail
*
* @param pos ...
* @return void
*/
void Airplane::addContactPoint(float* pos)
@ -638,8 +594,9 @@ void Airplane::compile()
t->handle = body->addMass(0, t->pos);
totalFuel += t->cap;
}
_cruiseWeight = _emptyWeight + totalFuel*_cruiseFuel;
_approachWeight = _emptyWeight + totalFuel*_approachFuel;
_cruiseConfig.weight = _emptyWeight + totalFuel*_cruiseConfig.fuel;
_approachConfig.weight = _emptyWeight + totalFuel*_approachConfig.fuel;
body->recalc();
@ -723,12 +680,12 @@ void Airplane::solveGear()
// The force at max compression should be sufficient to stop a
// plane moving downwards at 2x the approach descent rate. Assume
// a 3 degree approach.
float descentRate = 2.0f*_approachSpeed/19.1f;
float descentRate = 2.0f*_approachConfig.speed/19.1f;
// Spread the kinetic energy according to the gear weights. This
// will results in an equal compression fraction (not distance) of
// each gear.
float energy = 0.5f*_approachWeight*descentRate*descentRate;
float energy = 0.5f*_approachConfig.weight*descentRate*descentRate;
for(i=0; i<_gears.size(); i++) {
GearRec* gr = (GearRec*)_gears.get(i);
@ -743,7 +700,7 @@ void Airplane::solveGear()
gr->gear->setSpring(k * gr->gear->getSpring());
// Critically damped (too damped, too!)
gr->gear->setDamping(2*Math::sqrt(k*_approachWeight*gr->wgt)
gr->gear->setDamping(2*Math::sqrt(k*_approachConfig.weight*gr->wgt)
* gr->gear->getDamping());
}
}
@ -789,7 +746,7 @@ void Airplane::setupWeights(bool isApproach)
}
/// load values for controls as defined in cruise configuration
void Airplane::loadControls(Vector& controls)
void Airplane::loadControls(const Vector& controls)
{
_controls.reset();
for(int i=0; i < controls.size(); i++) {
@ -800,81 +757,43 @@ void Airplane::loadControls(Vector& controls)
}
/// Helper for solve()
void Airplane::runCruise()
void Airplane::runConfig(Config &cfg)
{
setupState(_cruiseAoA, _cruiseSpeed,_cruiseGlideAngle, &_cruiseState);
_model.setState(&_cruiseState);
_model.setAir(_cruiseP, _cruiseT,
Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
// The control configuration
loadControls(_cruiseControls);
// aoa is consider to be given for approach so we calculate orientation
// only once in setApproach()
if (!cfg.isApproach) {
cfg.state.setupOrientationFromAoa(cfg.aoa);
}
cfg.state.setupSpeedAndPosition(cfg.speed, cfg.glideAngle);
_model.setState(&cfg.state);
_model.setStandardAtmosphere(cfg.altitude);
loadControls(cfg.controls);
// The local wind
float wind[3];
Math::mul3(-1, _cruiseState.v, wind);
Math::vmul33(_cruiseState.orient, wind, wind);
Math::mul3(-1, cfg.state.v, wind);
cfg.state.globalToLocal(wind, wind);
setFuelFraction(_cruiseFuel);
setupWeights(false);
setFuelFraction(cfg.fuel);
setupWeights(cfg.isApproach);
// Set up the thruster parameters and iterate until the thrust
// stabilizes.
for(int i=0; i<_thrusters.size(); i++) {
Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
t->setWind(wind);
t->setAir(_cruiseP, _cruiseT,
Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
t->setStandardAtmosphere(cfg.altitude);
}
stabilizeThrust();
stabilizeThrust();
updateGearState();
// Precompute thrust in the model, and calculate aerodynamic forces
_model.getBody()->recalc();
_model.getBody()->reset();
_model.initIteration();
_model.calcForces(&_cruiseState);
_model.calcForces(&cfg.state);
}
/// Helper for solve()
void Airplane::runApproach()
{
setupState(_approachAoA, _approachSpeed,_approachGlideAngle, &_approachState);
_model.setState(&_approachState);
_model.setAir(_approachP, _approachT,
Atmosphere::calcStdDensity(_approachP, _approachT));
// The control configuration
loadControls(_approachControls);
// The local wind
float wind[3];
Math::mul3(-1, _approachState.v, wind);
Math::vmul33(_approachState.orient, wind, wind);
setFuelFraction(_approachFuel);
setupWeights(true);
// Run the thrusters until they get to a stable setting. FIXME:
// this is lots of wasted work.
for(int i=0; i<_thrusters.size(); i++) {
Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
t->setWind(wind);
t->setAir(_approachP, _approachT,
Atmosphere::calcStdDensity(_approachP, _approachT));
}
stabilizeThrust();
updateGearState();
// Precompute thrust in the model, and calculate aerodynamic forces
_model.getBody()->recalc();
_model.getBody()->reset();
_model.initIteration();
_model.calcForces(&_approachState);
}
/// Used only in Airplane::solve() and solveHelicopter(), not at runtime
void Airplane::applyDragFactor(float factor)
{
@ -960,51 +879,51 @@ void Airplane::solve()
}
// Run an iteration at cruise, and extract the needed numbers:
runCruise();
runConfig(_cruiseConfig);
_model.getThrust(tmp);
float thrust = tmp[0] + _cruiseWeight * Math::sin(_cruiseGlideAngle) * 9.81;
float thrust = tmp[0] + _cruiseConfig.weight * Math::sin(_cruiseConfig.glideAngle) * 9.81;
_model.getBody()->getAccel(tmp);
Math::tmul33(_cruiseState.orient, tmp, tmp);
float xforce = _cruiseWeight * tmp[0];
float clift0 = _cruiseWeight * tmp[2];
_cruiseConfig.state.localToGlobal(tmp, tmp);
float xforce = _cruiseConfig.weight * tmp[0];
float clift0 = _cruiseConfig.weight * tmp[2];
_model.getBody()->getAngularAccel(tmp);
Math::tmul33(_cruiseState.orient, tmp, tmp);
_cruiseConfig.state.localToGlobal(tmp, tmp);
float pitch0 = tmp[1];
// Run an approach iteration, and do likewise
runApproach();
runConfig(_approachConfig);
_model.getBody()->getAngularAccel(tmp);
Math::tmul33(_approachState.orient, tmp, tmp);
_approachConfig.state.localToGlobal(tmp, tmp);
double apitch0 = tmp[1];
_model.getBody()->getAccel(tmp);
Math::tmul33(_approachState.orient, tmp, tmp);
float alift = _approachWeight * tmp[2];
_approachConfig.state.localToGlobal(tmp, tmp);
float alift = _approachConfig.weight * tmp[2];
// Modify the cruise AoA a bit to get a derivative
_cruiseAoA += ARCMIN;
runCruise();
_cruiseAoA -= ARCMIN;
_cruiseConfig.aoa += ARCMIN;
runConfig(_cruiseConfig);
_cruiseConfig.aoa -= ARCMIN;
_model.getBody()->getAccel(tmp);
Math::tmul33(_cruiseState.orient, tmp, tmp);
float clift1 = _cruiseWeight * tmp[2];
_cruiseConfig.state.localToGlobal(tmp, tmp);
float clift1 = _cruiseConfig.weight * tmp[2];
// Do the same with the tail incidence
_tail->setIncidence(_tailIncidence + ARCMIN);
runCruise();
runConfig(_cruiseConfig);
_tail->setIncidence(_tailIncidence);
_model.getBody()->getAngularAccel(tmp);
Math::tmul33(_cruiseState.orient, tmp, tmp);
_cruiseConfig.state.localToGlobal(tmp, tmp);
float pitch1 = tmp[1];
// Now calculate:
float awgt = 9.8f * _approachWeight;
float awgt = 9.8f * _approachConfig.weight;
float dragFactor = thrust / (thrust-xforce);
float liftFactor = awgt / (awgt+alift);
@ -1021,11 +940,11 @@ void Airplane::solve()
// variable).
const float ELEVDIDDLE = 0.001f;
_approachElevator.val += ELEVDIDDLE;
runApproach();
runConfig(_approachConfig);
_approachElevator.val -= ELEVDIDDLE;
_model.getBody()->getAngularAccel(tmp);
Math::tmul33(_approachState.orient, tmp, tmp);
_approachConfig.state.localToGlobal(tmp, tmp);
double apitch1 = tmp[1];
float elevDelta = -apitch0 * (ELEVDIDDLE/(apitch1-apitch0));
@ -1044,14 +963,14 @@ void Airplane::solve()
}
// OK, now we can adjust the minor variables:
_cruiseAoA += SOLVE_TWEAK*aoaDelta;
_cruiseConfig.aoa += SOLVE_TWEAK*aoaDelta;
_tailIncidence += SOLVE_TWEAK*tailDelta;
_cruiseAoA = Math::clamp(_cruiseAoA, -0.175f, 0.175f);
_cruiseConfig.aoa = Math::clamp(_cruiseConfig.aoa, -0.175f, 0.175f);
_tailIncidence = Math::clamp(_tailIncidence, -0.175f, 0.175f);
if(abs(xforce/_cruiseWeight) < STHRESH*0.0001 &&
abs(alift/_approachWeight) < STHRESH*0.0001 &&
if(abs(xforce/_cruiseConfig.weight) < STHRESH*0.0001 &&
abs(alift/_approachConfig.weight) < STHRESH*0.0001 &&
abs(aoaDelta) < STHRESH*.000017 &&
abs(tailDelta) < STHRESH*.000017)
{
@ -1074,7 +993,7 @@ void Airplane::solve()
} else if(_liftRatio < 1e-04 || _liftRatio > 1e4) {
_failureMsg = "Lift ratio beyond reasonable bounds.";
return;
} else if(Math::abs(_cruiseAoA) >= .17453293) {
} else if(Math::abs(_cruiseConfig.aoa) >= .17453293) {
_failureMsg = "Cruise AoA > 10 degrees";
return;
} else if(Math::abs(_tailIncidence) >= .17453293) {
@ -1102,14 +1021,12 @@ void Airplane::solveHelicopter()
applyDragFactor(Math::pow(15.7/1000, 1/SOLVE_TWEAK));
applyLiftRatio(Math::pow(104, 1/SOLVE_TWEAK));
}
setupState(0,0,0, &_cruiseState);
_model.setState(&_cruiseState);
_cruiseConfig.state.setupState(0,0,0);
_model.setState(&_cruiseConfig.state);
setupWeights(true);
_controls.reset();
_model.getBody()->reset();
_model.setAir(_cruiseP, _cruiseT,
Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
_model.setStandardAtmosphere(_cruiseConfig.altitude);
}
float Airplane::getCGMAC()

126
src/FDM/YASim/Airplane.hpp
View file

@ -95,14 +95,13 @@ public:
int getSolutionIterations() const { return _solutionIterations; }
float getDragCoefficient() const { return _dragFactor; }
float getLiftRatio() const { return _liftRatio; }
float getCruiseAoA() const { return _cruiseAoA; }
float getCruiseAoA() const { return _cruiseConfig.aoa; }
float getTailIncidence() const { return _tailIncidence; }
float getApproachElevator() const { return _approachElevator.val; }
const char* getFailureMsg() const { return _failureMsg; }
static void setupState(const float aoa, const float speed, const float gla, yasim::State* s); // utility
void loadApproachControls() { loadControls(_approachControls); }
void loadCruiseControls() { loadControls(_cruiseControls); }
void loadApproachControls() { loadControls(_approachConfig.controls); }
void loadCruiseControls() { loadControls(_cruiseConfig.controls); }
float getCGHardLimitXMin() const { return _cgMin; } // get min x-coordinate for c.g (from main gear)
float getCGHardLimitXMax() const { return _cgMax; } // get max x-coordinate for c.g (from nose gear)
@ -114,21 +113,60 @@ public:
void setAutoBallast(bool allowed) { _autoBallast = allowed; }
private:
struct Tank { float pos[3]; float cap; float fill;
float density; int handle; };
struct Fuselage { float front[3], back[3], width, taper, mid, _cx, _cy, _cz, _idrag;
Vector surfs; };
struct GearRec { Gear* gear; Surface* surf; float wgt; };
struct ThrustRec { Thruster* thruster;
int handle; float cg[3]; float mass; };
struct Control { int control; float val; };
struct WeightRec { int handle; Surface* surf; };
struct SolveWeight { bool approach; int idx; float wgt; };
struct ContactRec { Gear* gear; float p[3]; };
struct Tank {
float pos[3];
float cap, fill, density;
int handle;
};
struct Fuselage {
float front[3], back[3];
float width, taper, mid, _cx, _cy, _cz, _idrag;
Vector surfs;
};
struct GearRec {
Gear* gear;
Surface* surf;
float wgt;
};
struct ThrustRec {
Thruster* thruster;
int handle;
float cg[3];
float mass;
};
struct Control {
int control;
float val;
};
struct WeightRec {
int handle;
Surface* surf;
};
struct SolveWeight {
bool approach;
int idx;
float wgt;
};
struct ContactRec {
Gear* gear;
float p[3];
};
struct Config {
bool isApproach {false};
float speed {0};
float fuel {0};
float glideAngle {0};
float aoa {0};
float altitude;
float weight;
State state;
Vector controls;
};
Config _cruiseConfig;
Config _approachConfig;
void loadControls(Vector &controls);
void runCruise();
void runApproach();
void loadControls(const Vector& controls);
void runConfig(Config &cfg);
void solveGear();
void solve();
void solveHelicopter();
@ -148,17 +186,17 @@ private:
Model _model;
ControlMap _controls;
float _emptyWeight;
float _pilotPos[3];
float _emptyWeight {0};
float _pilotPos[3] {0, 0, 0};
Wing* _wing;
Wing* _tail;
Wing* _wing {nullptr};
Wing* _tail {nullptr};
Vector _fuselages;
Vector _vstabs;
Vector _tanks;
Vector _thrusters;
float _ballast;
float _ballast {0};
Vector _gears;
Vector _contacts; // non-gear ground contact points
@ -168,39 +206,19 @@ private:
Vector _solveWeights;
Vector _cruiseControls;
State _cruiseState;
float _cruiseP;
float _cruiseT;
float _cruiseSpeed;
float _cruiseWeight;
float _cruiseFuel;
float _cruiseGlideAngle;
Vector _approachControls;
State _approachState;
float _approachP;
float _approachT;
float _approachSpeed;
float _approachAoA;
float _approachWeight;
float _approachFuel;
float _approachGlideAngle;
int _solutionIterations;
float _dragFactor;
float _liftRatio;
float _cruiseAoA;
float _tailIncidence;
int _solutionIterations {0};
float _dragFactor {1};
float _liftRatio {1};
float _tailIncidence {0};
Control _approachElevator;
const char* _failureMsg;
const char* _failureMsg {0};
float _cgMax = -1e6; // hard limits for cg from gear position
float _cgMin = 1e6; // hard limits for cg from gear position
float _cgDesiredMax = 0.3f; // desired cg max in %MAC from config
float _cgDesiredMin = 0.25f; // desired cg min in %MAC from config
float _cgDesiredFront; // calculated desired cg x max
float _cgDesiredAft; // calculated desired cg x min
float _cgMax {-1e6}; // hard limits for cg from gear position
float _cgMin {1e6}; // hard limits for cg from gear position
float _cgDesiredMax {0.3f}; // desired cg max in %MAC from config
float _cgDesiredMin {0.25f}; // desired cg min in %MAC from config
float _cgDesiredFront {0}; // calculated desired cg x max
float _cgDesiredAft {0}; // calculated desired cg x min
bool _autoBallast = false;
};

81
src/FDM/YASim/Atmosphere.cpp
View file

@ -1,7 +1,8 @@
#include <string>
#include "Math.hpp"
#include "Atmosphere.hpp"
namespace yasim {
namespace yasim {
// Copied from McCormick, who got it from "The ARDC Model Atmosphere"
// Note that there's an error in the text in the first entry,
// McCormick lists 299.16/101325/1.22500, but those don't agree with
@ -9,7 +10,7 @@ namespace yasim {
// pretty hot for a "standard" atmosphere.
// Numbers above 19000 meters calculated from src/Environment/environment.cxx
// meters kelvin Pa kg/m^3
float Atmosphere::data[][4] = {{ -900.0f, 293.91f, 111679.0f, 1.32353f },
float Atmosphere::data[][Atmosphere::numColumns] = {{ -900.0f, 293.91f, 111679.0f, 1.32353f },
{ 0.0f, 288.11f, 101325.0f, 1.22500f },
{ 900.0f, 282.31f, 90971.0f, 1.12260f },
{ 1800.0f, 276.46f, 81494.0f, 1.02690f },
@ -48,28 +49,34 @@ float Atmosphere::data[][4] = {{ -900.0f, 293.91f, 111679.0f, 1.32353f },
// Universal gas constant for air, in SI units. P = R * rho * T.
// P in pascals (N/m^2), rho is kg/m^3, T in kelvin.
const float R = 287.1f;
const float R = 287.058f;
// Specific heat ratio for air, at "low" temperatures.
const float GAMMA = 1.4f;
void Atmosphere::setStandard(float altitude)
{
_density = getStdDensity(altitude);
_pressure = getStdPressure(altitude);
_temperature = getStdTemperature(altitude);
}
float Atmosphere::getStdTemperature(float alt)
{
return getRecord(alt, 1);
return getRecord(alt, TEMPERATURE);
}
float Atmosphere::getStdPressure(float alt)
{
return getRecord(alt, 2);
return getRecord(alt, PRESSURE);
}
float Atmosphere::getStdDensity(float alt)
{
return getRecord(alt, 3);
return getRecord(alt, DENSITY);
}
float Atmosphere::calcVEAS(float spd,
float pressure, float temp, float density)
float Atmosphere::calcVEAS(float spd, float pressure, float temp, float density)
{
static float rho0 = getStdDensity(0);
float densityRatio = density / rho0;
@ -122,6 +129,11 @@ float Atmosphere::spdFromMach(float mach, float temp)
return mach * Math::sqrt(GAMMA * R * temp);
}
float Atmosphere::spdFromMach(float mach)
{
return spdFromMach(mach, _temperature);
}
float Atmosphere::spdFromVCAS(float vcas, float pressure, float temp)
{
// FIXME: does not account for supersonic
@ -136,6 +148,11 @@ float Atmosphere::spdFromVCAS(float vcas, float pressure, float temp)
return vtas;
}
float Atmosphere::spdFromVCAS(float vcas)
{
return spdFromVCAS(vcas, _pressure, _temperature);
}
void Atmosphere::calcStaticAir(float p0, float t0, float d0, float v,
float* pOut, float* tOut, float* dOut)
{
@ -147,14 +164,24 @@ void Atmosphere::calcStaticAir(float p0, float t0, float d0, float v,
*pOut = (*dOut) * R * (*tOut);
}
float Atmosphere::getRecord(float alt, int recNum)
void Atmosphere::calcStaticAir(float v, float* pOut, float* tOut, float* dOut)
{
int hi = (sizeof(data) / (4*sizeof(float))) - 1;
return calcStaticAir(_pressure, _temperature, _density, v, pOut, tOut, dOut);
}
float Atmosphere::getRecord(float alt, Column recNum)
{
int hi = maxTableIndex();
int lo = 0;
// safety valve, clamp to the edges of the table
if(alt < data[0][0]) hi=1;
else if(alt > data[hi][0]) lo = hi-1;
if(alt < data[0][ALTITUDE]) {
hi = 1;
}
else if(alt > data[hi][ALTITUDE]) {
lo = hi-1;
}
// binary search
while(1) {
@ -165,10 +192,38 @@ float Atmosphere::getRecord(float alt, int recNum)
}
// interpolate
float frac = (alt - data[lo][0])/(data[hi][0] - data[lo][0]);
float frac = (alt - data[lo][ALTITUDE])/(data[hi][ALTITUDE] - data[lo][ALTITUDE]);
float a = data[lo][recNum];
float b = data[hi][recNum];
return a + frac * (b-a);
}
int Atmosphere::maxTableIndex() {
return (sizeof(data) / (numColumns * sizeof(float))) - 1;
}
bool Atmosphere::test() {
bool passed = true;
int rows = maxTableIndex() + 1;
const float maxDeviation = 0.0002f;
fprintf(stderr, "Atmosphere::test()\n");
fprintf(stderr, "Columns = %d\n", numColumns);
fprintf(stderr, "Rows = %d\n", rows);
for (int alt = 0; alt < maxTableIndex(); alt++) {
float density = calcStdDensity(data[alt][PRESSURE], data[alt][TEMPERATURE]);
float delta = data[alt][DENSITY] - density;
fprintf(stderr, "%d : %f \n", alt, delta);
if (Math::abs(delta) > maxDeviation) {
passed = false;
fprintf(stderr,"FAIL: Deviation above limit of %1.6f\n", maxDeviation);
}
}
if (passed) {
fprintf(stderr,"Deviation below %1.6f for all rows.\n", maxDeviation);
}
return passed;
}
}; // namespace yasim

33
src/FDM/YASim/Atmosphere.hpp
View file

@ -3,8 +3,28 @@
namespace yasim {
//constexpr int Atmosphere::numColumns {4};
class Atmosphere {
enum Column {
ALTITUDE,
TEMPERATURE,
PRESSURE,
DENSITY
};
static const int numColumns {4};
public:
void setTemperature(float t) { _temperature = t; }
void setPressure(float p) { _pressure = p; }
void setDensity(float d) { _density = d; }
//set temperature, pressure and density to standard values for given altitude
void setStandard(float altitude);
float getTemperature() const { return _temperature; }
float getPressure() const { return _pressure; }
float getDensity() const { return _density; }
static float getStdTemperature(float alt);
static float getStdPressure(float alt);
static float getStdDensity(float alt);
@ -16,6 +36,8 @@ public:
static float spdFromMach(float mach, float temp);
static float spdFromVCAS(float vcas, float pressure, float temp);
float spdFromMach(float mach);
float spdFromVCAS(float vcas);
// Given ambient ("0") pressure/density/temperature values,
// calculate the properties of static air (air accelerated to the
@ -23,10 +45,17 @@ public:
// compressibility, but not shock effects.
static void calcStaticAir(float p0, float t0, float d0, float v,
float* pOut, float* tOut, float* dOut);
void calcStaticAir(float v, float* pOut, float* tOut, float* dOut);
static bool test();
private:
static float getRecord(float alt, int idx);
static float data[][4];
static float getRecord(float alt, Atmosphere::Column recNum);
static float data[][numColumns];
static int maxTableIndex();
float _temperature = 288.11f;
float _pressure = 101325.0f;
float _density = 1.22500f;
};
}; // namespace yasim

66
src/FDM/YASim/BodyEnvironment.hpp
View file

@ -13,45 +13,41 @@ namespace yasim {
// You can get local->global transformations by calling Math::tmul33()
// and using the same matrix.
struct State {
double pos[3]; // position
float orient[9]; // global->local xform matrix
float v[3]; // velocity
float rot[3]; // rotational velocity
float acc[3]; // acceleration
float racc[3]; // rotational acceleration
double pos[3] {0, 0, 0}; // position
float orient[9] {1,0,0, 0,1,0, 0,0,1}; // global->local xform matrix
float v[3] {0, 0, 0}; // velocity
float rot[3] {0, 0, 0}; // rotational velocity
float acc[3] {0, 0, 0}; // acceleration
float racc[3] {0, 0, 0}; // rotational acceleration
// Simple initialization
State() {
int i;
for(i=0; i<3; i++) {
pos[i] = v[i] = rot[i] = acc[i] = racc[i] = 0;
int j;
for(j=0; j<3; j++)
orient[3*i+j] = i==j ? 1.0f : 0.0f;
}
}
void posLocalToGlobal(float* lpos, double *gpos) {
void posLocalToGlobal(const float* lpos, double *gpos) const {
float tmp[3];
Math::tmul33(orient, lpos, tmp);
gpos[0] = tmp[0] + pos[0];
gpos[1] = tmp[1] + pos[1];
gpos[2] = tmp[2] + pos[2];
}
void posGlobalToLocal(double* gpos, float *lpos) {
void posGlobalToLocal(const double* gpos, float *lpos) const {
lpos[0] = (float)(gpos[0] - pos[0]);
lpos[1] = (float)(gpos[1] - pos[1]);
lpos[2] = (float)(gpos[2] - pos[2]);
Math::vmul33(orient, lpos, lpos);
}
void velLocalToGlobal(float* lvel, float *gvel) {
Math::tmul33(orient, lvel, gvel);
}
void velGlobalToLocal(float* gvel, float *lvel) {
Math::vmul33(orient, gvel, lvel);
void localToGlobal(const float* local, float *global) const {
Math::tmul33(orient, local, global);
}
void planeGlobalToLocal(double* gplane, float *lplane) {
void globalToLocal(const float* global, float *local) const {
Math::vmul33(orient, global, local);
}
void planeGlobalToLocal(const double* gplane, float *lplane) const {
// First the normal vector transformed to local coordinates.
lplane[0] = (float)-gplane[0];
lplane[1] = (float)-gplane[1];
@ -63,6 +59,34 @@ struct State {
+ pos[2]*gplane[2] - gplane[3]);
}
// used by Airplane::runCruise, runApproach, solveHelicopter and in yasim-test
void setupOrientationFromAoa(float aoa)
{
float cosAoA = Math::cos(aoa);
float sinAoA = Math::sin(aoa);
orient[0] = cosAoA; orient[1] = 0; orient[2] = sinAoA;
orient[3] = 0; orient[4] = 1; orient[5] = 0;
orient[6] = -sinAoA; orient[7] = 0; orient[8] = cosAoA;
}
void setupSpeedAndPosition(float speed, float gla)
{
// FIXME check axis, guess sin should go to 2 instead of 1?
v[0] = speed*Math::cos(gla);
v[1] = -speed*Math::sin(gla);
v[2] = 0;
for(int i=0; i<3; i++) {
pos[i] = rot[i] = acc[i] = racc[i] = 0;
}
pos[2] = 1;
}
void setupState(float aoa, float speed, float gla) {
setupOrientationFromAoa(aoa);
setupSpeedAndPosition(speed, gla);
}
};
//

2
src/FDM/YASim/CMakeLists.txt
View file

@ -22,7 +22,6 @@ set(COMMON
Rotorpart.cpp
SimpleJet.cpp
Surface.cpp
Thruster.cpp
TurbineEngine.cpp
Turbulence.cpp
Wing.cpp
@ -40,6 +39,7 @@ flightgear_component(YASim "${SOURCES}")
if(ENABLE_TESTS)
add_executable(yasim yasim-test.cpp ${COMMON})
add_executable(yasim-proptest proptest.cpp ${COMMON})
add_executable(yasim-atmotest yasim-atmotest.cpp Atmosphere.cpp )
target_link_libraries(yasim SimGearCore)
target_link_libraries(yasim-proptest SimGearCore)

19
src/FDM/YASim/FGFDM.cpp
View file

@ -27,25 +27,6 @@
namespace yasim {
// Some conversion factors
static const float KTS2MPS = 0.514444444444;
static const float KMH2MPS = 1/3.6;
static const float FT2M = 0.3048;
static const float DEG2RAD = 0.0174532925199;
static const float RPM2RAD = 0.10471975512;
static const float LBS2N = 4.44822;
static const float LBS2KG = 0.45359237;
static const float KG2LBS = 2.2046225;
static const float CM2GALS = 264.172037284;
static const float HP2W = 745.700;
static const float INHG2PA = 3386.389;
static const float K2DEGF = 1.8;
static const float K2DEGFOFFSET = -459.4;
static const float CIN2CM = 1.6387064e-5;
static const float YASIM_PI = 3.14159265358979323846;
static const float NM2FTLB = (1/(LBS2N*FT2M));
// Stubs, so that this can be compiled without the FlightGear
// binary. What's the best way to handle this?

1
src/FDM/YASim/FGFDM.hpp
View file

@ -4,6 +4,7 @@
#include <simgear/xml/easyxml.hxx>
#include <simgear/props/props.hxx>
#include "yasim-common.hpp"
#include "Airplane.hpp"
#include "Vector.hpp"

2
src/FDM/YASim/Gear.cpp
View file

@ -165,7 +165,7 @@ void Gear::calcForce(RigidBody* body, State *s, float* v, float* rot)
// The velocity of the contact patch transformed to local coordinates.
float glvel[3];
s->velGlobalToLocal(_global_vel, glvel);
s->globalToLocal(_global_vel, glvel);
// First off, make sure that the gear "tip" is below the ground.
// If it's not, there's no force.

6
src/FDM/YASim/Hitch.cpp
View file

@ -413,7 +413,7 @@ void Hitch::calcForce(Ground *g_cb, RigidBody* body, State* s)
//With this trick, both player in aerotow get the same length
Math::unit3(delta,deltaN);
float lvel[3];
s->velGlobalToLocal(s->v,lvel);
s->globalToLocal(s->v,lvel);
_speed_in_tow_direction=Math::dot3(lvel,deltaN);
if (_towEndIsConnectedToProperty && _nodeIsMultiplayer)
{
@ -421,7 +421,7 @@ void Hitch::calcForce(Ground *g_cb, RigidBody* body, State* s)
_timeLagCorrectedDist=_dist+mp_delta_dist_due_to_time_lag;
if(_forceIsCalculatedByMaster && !_open)
{
s->velGlobalToLocal(_mp_force,_force);
s->globalToLocal(_mp_force,_force);
return;
}
}
@ -528,7 +528,7 @@ void Hitch::calcForce(Ground *g_cb, RigidBody* body, State* s)
//the same for the tow end:
Math::mul3(grav_frac_tow_end*grav_force,ground,grav_force_v);
Math::add3(grav_force_v,_towEndForce,_towEndForce);
s->velLocalToGlobal(_towEndForce,_towEndForce);
s->localToGlobal(_towEndForce,_towEndForce);
if(_forceMagnitude>=_towBrakeForce)
{

4
src/FDM/YASim/Hook.cpp
View file

@ -216,8 +216,8 @@ void Hook::calcForce(Ground* g_cb, RigidBody* body, State* s, float* lv, float*
float wire_lpos[2][3];
s->posGlobalToLocal(dpos[0], wire_lpos[0]);
s->posGlobalToLocal(dpos[1], wire_lpos[1]);
s->velGlobalToLocal(wire_vel[0], wire_vel[0]);
s->velGlobalToLocal(wire_vel[1], wire_vel[1]);
s->globalToLocal(wire_vel[0], wire_vel[0]);
s->globalToLocal(wire_vel[1], wire_vel[1]);
// Compute the velocity of the hooks mount point in the local frame.
float mount_vel[3];

37
src/FDM/YASim/Integrator.cpp
View file

@ -2,43 +2,6 @@
#include "Integrator.hpp"
namespace yasim {
void Integrator::setBody(RigidBody* body)
{
_body = body;
}
void Integrator::setEnvironment(BodyEnvironment* env)
{
_env = env;
}
void Integrator::setInterval(float dt)
{
_dt = dt;
}
float Integrator::getInterval()
{
return _dt;
}
void Integrator::setState(State* s)
{
_s = *s;
}
State* Integrator::getState()
{
return &_s;
}
// Transforms a "local" vector to a "global" vector (not coordinate!)
// using the specified orientation.
void Integrator::l2gVector(float* orient, float* v, float* out)
{
Math::tmul33(orient, v, out);
}
// Updates a position vector for a body c.g. motion with velocity v
// over time dt, from orientation o0 to o1. Because the position
// references the local coordinate origin, but the velocity is that of

22
src/FDM/YASim/Integrator.hpp
View file

@ -17,24 +17,24 @@ class Integrator
{
public:
// Sets the RigidBody that will be integrated.
void setBody(RigidBody* body);
void setBody(RigidBody* body) { _body = body; }
// Sets the BodyEnvironment object used to calculate the second
// derivatives.
void setEnvironment(BodyEnvironment* env);
void setEnvironment(BodyEnvironment* env) { _env = env; }
// Sets the time interval between steps. Units can be anything,
// but they must match the other values (if you specify speed in
// m/s, then angular acceleration had better be in 1/s^2 and the
// time interval should be in seconds, etc...)
void setInterval(float dt);
float getInterval();
void setInterval(float dt) { _dt = dt; }
float getInterval() const { return _dt; }
// The current state, i.e. initial conditions for the next
// integration iteration. Note that the acceleration parameters
// in the State object are ignored.
State* getState();
void setState(State* s);
State* getState() { return &_s; }
void setState(State* s) { _s = *s; }
// Do a 4th order Runge-Kutta integration over one time interval.
// This is the top level of the simulation.
@ -43,9 +43,13 @@ public:
private:
void orthonormalize(float* m);
void rotMatrix(float* r, float dt, float* out);
void l2gVector(float* orient, float* v, float* out);
void extrapolatePosition(double* pos, float* v, float dt,
float* o1, float* o2);
// Transforms a "local" vector to a "global" vector (not coordinate!)
// using the specified orientation.
void l2gVector(const float* orient, const float* v, float* out) {
Math::tmul33(orient, v, out); }
void extrapolatePosition(double* pos, float* v, float dt, float* o1, float* o2);
BodyEnvironment* _env;
RigidBody* _body;

65
src/FDM/YASim/Jet.cpp
View file

@ -59,21 +59,6 @@ void Jet::setMaxThrust(float thrust, float afterburner)
}
}
void Jet::setVMax(float spd)
{
_vMax = spd;
}
void Jet::setTSFC(float tsfc)
{
_tsfc = tsfc;
}
void Jet::setATSFC(float atsfc)
{
_atsfc = atsfc;
}
void Jet::setRPMs(float idleN1, float maxN1, float idleN2, float maxN2)
{
_n1Min = idleN1;
@ -82,16 +67,6 @@ void Jet::setRPMs(float idleN1, float maxN1, float idleN2, float maxN2)
_n2Max = maxN2;
}
void Jet::setEGT(float takeoffEGT)
{
_egt0 = takeoffEGT;
}
void Jet::setEPR(float takeoffEPR)
{
_epr0 = takeoffEPR;
}
void Jet::setSpooling(float time)
{
// 2.3 = -ln(0.1), i.e. x=2.3 is the 90% point we're defining
@ -100,11 +75,6 @@ void Jet::setSpooling(float time)
_decay = 1.5f * 2.3f / time;
}
void Jet::setVectorAngle(float angle)
{
_maxRot = angle;
}
void Jet::setReheat(float reheat)
{
_reheat = Math::clamp(reheat, 0, 1);
@ -112,9 +82,7 @@ void Jet::setReheat(float reheat)
void Jet::setRotation(float rot)
{
if(rot < 0) rot = 0;
if(rot > 1) rot = 1;
_rotControl = rot;
_rotControl = Math::clamp(rot, 0, 1);
}
float Jet::getN1()
@ -127,16 +95,12 @@ float Jet::getN2()
return _n2 * _tempCorrect;
}
float Jet::getEPR()
{
return _epr;
}
float Jet::getEGT()
{
// Exactly zero means "off" -- return the ambient temperature
if(_egt == 0) return _temp;
if(_egt == 0) {
return _atmo.getTemperature();
}
return _egt * _tempCorrect * _tempCorrect;
}
@ -155,8 +119,7 @@ void Jet::integrate(float dt)
float spd = -Math::dot3(_wind, _dir);
float statT, statP, statD;
Atmosphere::calcStaticAir(_pressure, _temp, _rho, spd,
&statP, &statT, &statD);
_atmo.calcStaticAir(spd, &statP, &statT, &statD);
_pressureCorrect = statP/P0;
_tempCorrect = Math::sqrt(statT/T0);
@ -175,8 +138,8 @@ void Jet::integrate(float dt)
// Now get a "beta" (i.e. EPR - 1) value. The output values are
// expressed as functions of beta.
float ibeta0 = 1/(_epr0 - 1);
float betaTarget = (_epr0 - 1) * (setThrust/_maxThrust) * (P0/_pressure)
* (_temp/statT);
float betaTarget = (_epr0 - 1) * (setThrust/_maxThrust) * (P0/_atmo.getPressure())
* (_atmo.getTemperature()/statT);
float n2Target = _n2Min + (betaTarget*ibeta0) * (_n2Max - _n2Min);
// Note that this "first" beta value is used to compute a target
@ -191,7 +154,7 @@ void Jet::integrate(float dt)
// The actual thrust produced is keyed to the N1 speed. Add the
// afterburners in at the end.
float betaN1 = (_epr0-1) * (_n1 - _n1Min) / (_n1Max - _n1Min);
_thrust = _maxThrust * betaN1/((_epr0-1)*(P0/_pressure)*(_temp/statT));
_thrust = _maxThrust * betaN1/((_epr0-1)*(P0/_atmo.getPressure())*(_atmo.getTemperature()/statT));
_thrust *= 1 + _reheat*(_abFactor-1);
// Finally, calculate the output variables. Use a 80/20 mix of
@ -214,16 +177,6 @@ void Jet::integrate(float dt)
if(_reverseThrust) _thrust *= -_reverseEff;
}
bool Jet::isRunning()
{
return _running;
}
bool Jet::isCranking()
{
return _cranking;
}
void Jet::getThrust(float* out)
{
Math::mul3(_thrust, _dir, out);
@ -240,13 +193,11 @@ void Jet::getThrust(float* out)
void Jet::getTorque(float* out)
{
out[0] = out[1] = out[2] = 0;
return;
}
void Jet::getGyro(float* out)
{
out[0] = out[1] = out[2] = 0;
return;
}
}; // namespace yasim

18
src/FDM/YASim/Jet.hpp
View file

@ -12,13 +12,13 @@ public:
virtual Jet* getJet() { return this; }
void setMaxThrust(float thrust, float afterburner=0);
void setVMax(float spd);
void setTSFC(float tsfc);
void setATSFC(float atsfc);
void setVMax(float spd) { _vMax = spd; };
void setTSFC(float tsfc) { _tsfc = tsfc; };
void setATSFC(float atsfc) { _atsfc = atsfc; };
void setRPMs(float idleN1, float maxN1, float idleN2, float maxN2);
void setEGT(float takeoffEGT);
void setEPR(float takeoffEPR);
void setVectorAngle(float angle);
void setEGT(float takeoffEGT) { _egt0 = takeoffEGT; };
void setEPR(float takeoffEPR) { _epr0 = takeoffEPR; };
void setVectorAngle(float angle) { _maxRot = angle; };
// The time it takes the engine to reach 90% thrust from idle
void setSpooling(float time);
@ -37,15 +37,15 @@ public:
float getN1();
float getN2();
float getEPR();
float getEPR() const { return _epr; };
float getEGT();
// Normalized "performance" number. Used for fuzzy numbers in FGFDM
float getPerfNorm() { return (_n1 - _n1Min) / (_n1Max - _n1Min); }
// From Thruster:
virtual bool isRunning();
virtual bool isCranking();
virtual bool isRunning() { return _running; };
virtual bool isCranking() { return _cranking; };
virtual void getThrust(float* out);
virtual void getTorque(float* out);
virtual void getGyro(float* out);

4
src/FDM/YASim/Launchbar.cpp
View file

@ -374,8 +374,8 @@ void Launchbar::calcForce(Ground *g_cb, RigidBody* body, State* s, float* lv, fl
// Transform the velocities of the endpoints to the
// local coordinate sytem.
float lvel[2][3];
s->velGlobalToLocal(vel[0], lvel[0]);
s->velGlobalToLocal(vel[1], lvel[1]);
s->globalToLocal(vel[0], lvel[0]);
s->globalToLocal(vel[1], lvel[1]);
// Compute the position of the launchbar tip relative to the cat.
float tip_pos_on_cat = getPercentPosOnCat(llb_mount, 0.0, lend);

48
src/FDM/YASim/Model.cpp
View file

@ -50,28 +50,13 @@ void printState(State* s)
Model::Model()
{
int i;
for(i=0; i<3; i++) _wind[i] = 0;
_integrator.setBody(&_body);
_integrator.setEnvironment(this);
// Default value of 30 Hz
_integrator.setInterval(1.0f/30.0f);
_agl = 0;
_crashed = false;
_turb = 0;
_ground_cb = new Ground();
_hook = 0;
_launchbar = 0;
_wingSpan = 0;
_groundEffect = 0;
for(i=0; i<3; i++) _geRefPoint[i] = 0;
_global_ground[0] = 0; _global_ground[1] = 0; _global_ground[2] = 1;
_global_ground[3] = -100000;
_modelN = fgGetNode("/fdm/yasim/forces", true);
_fAeroXN = _modelN->getNode("f-x-drag", true);
_fAeroYN = _modelN->getNode("f-y-side", true);
@ -98,8 +83,7 @@ void Model::getThrust(float* out) const
{
float tmp[3];
out[0] = out[1] = out[2] = 0;
int i;
for(i=0; i<_thrusters.size(); i++) {
for(int i=0; i<_thrusters.size(); i++) {
Thruster* t = (Thruster*)_thrusters.get(i);
t->getThrust(tmp);
Math::add3(tmp, out, out);
@ -127,7 +111,7 @@ void Model::initIteration()
localWind(pos, _s, v, alt);
t->setWind(v);
t->setAir(_pressure, _temp, _rho);
t->setAtmosphere(_atmo);
t->integrate(_integrator.getInterval());
t->getTorque(v);
@ -153,8 +137,6 @@ void Model::initIteration()
Hitch* h = (Hitch*)_hitches.get(i);
h->integrate(_integrator.getInterval());
}
}
// This function initializes some variables for the rotor calculation
@ -186,7 +168,7 @@ void Model::iterate()
void Model::setState(State* s)
{
_integrator.setState(s);
_s = _integrator.getState();
_s = s;
}
@ -203,20 +185,6 @@ void Model::setGroundEffect(const float* pos, float span, float mul)
_groundEffect = mul;
}
void Model::setAir(float pressure, float temp, float density)
{
_pressure = pressure;
_temp = temp;
_rho = density;
}
void Model::setAirFromStandardAtmosphere(float altitude)
{
_pressure = Atmosphere::getStdPressure(altitude);
_temp = Atmosphere::getStdTemperature(altitude);
_rho = Atmosphere::getStdDensity(altitude);
}
void Model::updateGround(State* s)
{
float dummy[3];
@ -336,7 +304,7 @@ void Model::calcForces(State* s)
localWind(pos, s, vs, alt);
float force[3], torque[3];
sf->calcForce(vs, _rho, force, torque);
sf->calcForce(vs, _atmo.getDensity(), force, torque);
Math::add3(faero, force, faero);
_body.addForce(pos, force);
@ -349,7 +317,7 @@ void Model::calcForces(State* s)
float vs[3], pos[3];
r->getPosition(pos);
localWind(pos, s, vs, alt);
r->calcLiftFactor(vs, _rho,s);
r->calcLiftFactor(vs, _atmo.getDensity(), s);
float tq=0;
// total torque of rotor (scalar) for calculating new rotor rpm
@ -363,7 +331,7 @@ void Model::calcForces(State* s)
localWind(pos, s, vs, alt,true);
float force[3], torque[3];
rp->calcForce(vs, _rho, force, torque, &torque_scalar);
rp->calcForce(vs, _atmo.getDensity(), force, torque, &torque_scalar);
tq+=torque_scalar;
rp->getPositionForceAttac(pos);
@ -486,7 +454,7 @@ void Model::newState(State* s)
// Calculates the airflow direction at the given point and for the
// specified aircraft velocity.
void Model::localWind(const float* pos, yasim::State* s, float* out, float alt, bool is_rotor)
void Model::localWind(const float* pos, const yasim::State* s, float* out, float alt, bool is_rotor)
{
float tmp[3], lwind[3], lrot[3], lv[3];
@ -521,8 +489,6 @@ void Model::localWind(const float* pos, yasim::State* s, float* out, float alt,
_rotorgear.getDownWash(pos,lv,tmp);
Math::add3(out,tmp, out); // + downwash
}
}
}; // namespace yasim

36
src/FDM/YASim/Model.hpp
View file

@ -7,6 +7,7 @@
#include "Vector.hpp"
#include "Turbulence.hpp"
#include "Rotor.hpp"
#include "Atmosphere.hpp"
#include <simgear/props/props.hxx>
namespace yasim {
@ -68,8 +69,8 @@ public:
//
void setGroundEffect(const float* pos, float span, float mul);
void setWind(float* wind) { Math::set3(wind, _wind); }
void setAir(float pressure, float temp, float density);
void setAirFromStandardAtmosphere(float altitude);
void setAtmosphere(Atmosphere a) { _atmo = a; };
void setStandardAtmosphere(float altitude) { _atmo.setStandard(altitude); };
void updateGround(State* s);
@ -81,39 +82,38 @@ private:
void initRotorIteration();
void calcGearForce(Gear* g, float* v, float* rot, float* ground);
float gearFriction(float wgt, float v, Gear* g);
void localWind(const float* pos, State* s, float* out, float alt, bool is_rotor = false);
void localWind(const float* pos, const yasim::State* s, float* out, float alt, bool is_rotor = false);
Integrator _integrator;
RigidBody _body;
Turbulence* _turb;
Turbulence* _turb {nullptr};
Vector _thrusters;
Vector _surfaces;
Rotorgear _rotorgear;
Vector _gears;
Hook* _hook;
Launchbar* _launchbar;
Hook* _hook {nullptr};
Launchbar* _launchbar {nullptr};
Vector _hitches;
float _wingSpan;
float _groundEffect;
float _geRefPoint[3];
float _wingSpan {0};
float _groundEffect {0};
float _geRefPoint[3] {0,0,0};
Ground* _ground_cb;
double _global_ground[4];
float _pressure;
float _temp;
float _rho;
float _wind[3];
double _global_ground[4] {0,0,1, -1e5};
Atmosphere _atmo;
float _wind[3] {0,0,0};
// Accumulators for the total internal gyro and engine torque
float _gyro[3];
float _torque[3];
float _gyro[3] {0,0,0};
float _torque[3] {0,0,0};
State* _s;
bool _crashed;
float _agl;
bool _crashed {false};
float _agl {0};
SGPropertyNode_ptr _modelN;
SGPropertyNode_ptr _fAeroXN;
SGPropertyNode_ptr _fAeroYN;

10
src/FDM/YASim/PropEngine.cpp
View file

@ -123,8 +123,8 @@ void PropEngine::stabilize()
// If we cannot manage this in 100 iterations, give up.
for (int n = 0; n < 100; n++) {
float ptau, thrust;
_prop->calc(_rho, speed, _omega * _gearRatio, &thrust, &ptau);
_eng->calc(_pressure, _temp, _omega);
_prop->calc(_atmo.getDensity(), speed, _omega * _gearRatio, &thrust, &ptau);
_eng->calc(_atmo.getPressure(), _atmo.getTemperature(), _omega);
_eng->stabilize();
// Do it again -- the turbo sets the target MP in the first
@ -133,7 +133,7 @@ void PropEngine::stabilize()
// it works without side effects (other than solver
// performance). In the future, the Engine objects should
// store state to allow them to do the work themselves.
_eng->calc(_pressure, _temp, _omega);
_eng->calc(_atmo.getPressure(), _atmo.getTemperature(), _omega);
// Compute torque as seen by the engine's end of the gearbox.
// The propeller will be moving more slowly (for gear ratios
@ -185,11 +185,11 @@ void PropEngine::integrate(float dt)
_eng->setMixture(_mixture);
_eng->setFuelState(_fuel);
_prop->calc(_rho, speed, _omega * _gearRatio, &thrust, &propTorque);
_prop->calc(_atmo.getDensity(), speed, _omega * _gearRatio, &thrust, &propTorque);
if(_omega == 0.0)
_omega = 0.001; // hack to get around reports of NaNs somewhere...
propTorque *= _gearRatio;
_eng->calc(_pressure, _temp, _omega);
_eng->calc(_atmo.getPressure(), _atmo.getTemperature(), _omega);
_eng->integrate(dt);
engTorque = _eng->getTorque();
_fuelFlow = _eng->getFuelFlow();

2
src/FDM/YASim/Rotor.cpp
View file

@ -829,7 +829,7 @@ void Rotor::calcLiftFactor(float* v, float rho, State *s)
//store the gravity direction
Glue::geodUp(s->pos, _grav_direction);
s->velGlobalToLocal(_grav_direction, _grav_direction);
s->globalToLocal(_grav_direction, _grav_direction);
}
void Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s)

25
src/FDM/YASim/Surface.cpp
View file

@ -7,35 +7,12 @@ int Surface::s_idGenerator = 0;
Surface::Surface( Version * version ) :
_version(version)
{
// create id for surface
_id = s_idGenerator++;
// Start in a "sane" mode, so unset stuff doesn't freak us out
_c0 = 1;
_cx = _cy = _cz = 1;
_cz0 = 0;
_peaks[0] = _peaks[1] = 1;
int i;
for(i=0; i<4; i++) {
_stalls[i] = 0;
_widths[i] = 0.01; // half a degree
}
_orient[0] = 1; _orient[1] = 0; _orient[2] = 0;
_orient[3] = 0; _orient[4] = 1; _orient[5] = 0;
_orient[6] = 0; _orient[7] = 0; _orient[8] = 1;
_chord = 0;
_incidence = 0;
_twist = 0;
_slatPos = _spoilerPos = _flapPos = 0;
_slatDrag = _spoilerDrag = _flapDrag = 1;
_flapLift = 0;
_flapEffectiveness = 1;
_slatAlpha = 0;
_spoilerLift = 1;
_inducedDrag = 1;
_stallAlpha = 0;
_alpha = 0;
_surfN = fgGetNode("/fdm/yasim/debug/surfaces", true);
if (_surfN != 0) {
_surfN = _surfN->getChild("surface", _id, true);

48
src/FDM/YASim/Surface.hpp
View file

@ -92,36 +92,36 @@ private:
float flapLift(float alpha);
float controlDrag(float lift, float drag);
float _chord; // X-axis size
float _c0; // total force coefficient
float _cx; // X-axis force coefficient
float _cy; // Y-axis force coefficient
float _cz; // Z-axis force coefficient
float _cz0; // Z-axis force offset
float _peaks[2]; // Stall peak coefficients (fwd, back)
float _stalls[4]; // Stall angles (fwd/back, pos/neg)
float _widths[4]; // Stall widths " "
float _chord {0}; // X-axis size
float _c0 {1}; // total force coefficient
float _cx {1}; // X-axis force coefficient
float _cy {1}; // Y-axis force coefficient
float _cz {1}; // Z-axis force coefficient
float _cz0 {0}; // Z-axis force offset
float _peaks[2] {1, 1}; // Stall peak coefficients (fwd, back)
float _stalls[4] {0, 0, 0, 0}; // Stall angles (fwd/back, pos/neg)
float _widths[4] {0.01, 0.01, 0.01, 0.01}; // Stall widths " "
float _pos[3]; // position in local coords
float _orient[9]; // local->surface orthonormal matrix
float _slatAlpha;
float _slatDrag;
float _flapLift;
float _flapDrag;
float _flapEffectiveness;
float _spoilerLift;
float _spoilerDrag;
float _slatAlpha {0};
float _slatDrag {1};
float _flapLift {0};
float _flapDrag {1};
float _flapEffectiveness {1};
float _spoilerLift {1};
float _spoilerDrag {1};
float _slatPos;
float _flapPos;
float _spoilerPos;
float _incidence;
float _twist;
float _inducedDrag;
float _slatPos {0};
float _flapPos {0};
float _spoilerPos {0};
float _incidence {0};
float _twist {0};
float _inducedDrag {1};
// used during calculations
float _stallAlpha;
float _alpha;
float _stallAlpha {0};
float _alpha {0};
Version * _version;
SGPropertyNode* _fxN;

72
src/FDM/YASim/Thruster.cpp
View file

@ -1,72 +0,0 @@
#include "Math.hpp"
#include "Thruster.hpp"
namespace yasim {
Thruster::Thruster()
{
_dir[0] = 1; _dir[1] = 0; _dir[2] = 0;
int i;
for(i=0; i<3; i++) _pos[i] = _wind[i] = 0;
_throttle = 0;
_mixture = 0;
_starter = false;
_pressure = _temp = _rho = 0;
}
Thruster::~Thruster()
{
}
void Thruster::getPosition(float* out)
{
int i;
for(i=0; i<3; i++) out[i] = _pos[i];
}
void Thruster::setPosition(float* pos)
{
int i;
for(i=0; i<3; i++) _pos[i] = pos[i];
}
void Thruster::getDirection(float* out)
{
int i;
for(i=0; i<3; i++) out[i] = _dir[i];
}
void Thruster::setDirection(float* dir)
{
Math::unit3(dir, _dir);
}
void Thruster::setThrottle(float throttle)
{
_throttle = Math::clamp(throttle, -1, 1);
}
void Thruster::setMixture(float mixture)
{
_mixture = Math::clamp(mixture, 0, 1);
}
void Thruster::setStarter(bool starter)
{
_starter = starter;
}
void Thruster::setWind(float* wind)
{
int i;
for(i=0; i<3; i++) _wind[i] = wind[i];
}
void Thruster::setAir(float pressure, float temp, float density)
{
_pressure = pressure;
_temp = temp;
_rho = density;
}
}; // namespace yasim

43
src/FDM/YASim/Thruster.hpp
View file

@ -1,6 +1,9 @@
#ifndef _THRUSTER_HPP
#define _THRUSTER_HPP
#include "Atmosphere.hpp"
#include "Math.hpp"
namespace yasim {
class Jet;
@ -10,8 +13,8 @@ class Engine;
class Thruster {
public:
Thruster();
virtual ~Thruster();
Thruster() {};
virtual ~Thruster() {};
// Typing info, these are the possible sub-type (or sub-parts)
// that a thruster might have. Any might return null. A little
@ -22,15 +25,15 @@ public:
virtual Engine* getEngine() { return 0; }
// Static data
void getPosition(float* out);
void setPosition(float* pos);
void getDirection(float* out);
void setDirection(float* dir);
void getPosition(float* out) const { Math::set3(_pos, out); }
void setPosition(const float* pos) { Math::set3(pos, _pos); }
void getDirection(float* out) const { Math::set3(_dir, out); }
void setDirection(const float* dir) { Math::unit3(dir, _dir); }
// Controls
void setThrottle(float throttle);
void setMixture(float mixture);
void setStarter(bool starter);
void setThrottle(float throttle) { _throttle = Math::clamp(throttle, -1, 1); }
void setMixture(float mixture) { _mixture = Math::clamp(mixture, 0, 1); }
void setStarter(bool starter) { _starter = starter; }
void setFuelState(bool hasFuel) { _fuel = hasFuel; }
// Dynamic output
@ -42,26 +45,24 @@ public:
virtual float getFuelFlow()=0; // in kg/s
// Runtime instructions
void setWind(float* wind);
void setAir(float pressure, float temp, float density);
void setWind(const float* wind) { Math::set3(wind, _wind); };
void setAtmosphere(Atmosphere a) { _atmo = a; };
void setStandardAtmosphere(float altitude) { _atmo.setStandard(altitude); };
virtual void init() {}
virtual void integrate(float dt)=0;
virtual void stabilize()=0;
protected:
float _pos[3];
float _dir[3];
float _throttle;
float _mixture;
bool _starter; // true=engaged, false=disengaged
float _pos[3] {0, 0, 0};
float _dir[3] {1, 0, 0};
float _throttle = 0;
float _mixture = 0;
bool _starter = false; // true=engaged, false=disengaged
bool _fuel; // true=available, false=out
float _wind[3];
float _pressure;
float _temp;
float _rho;
float _wind[3] {0, 0, 0};
Atmosphere _atmo;
};
}; // namespace yasim
#endif // _THRUSTER_HPP

10
src/FDM/YASim/Version.hpp
View file

@ -16,20 +16,20 @@ public:
} YASIM_VERSION;
void setVersion( const char * version );
int getVersion() { return _version; }
bool isVersion( YASIM_VERSION version );
bool isVersionOrNewer( YASIM_VERSION version );
int getVersion() const { return _version; }
bool isVersion( YASIM_VERSION version ) const;
bool isVersionOrNewer( YASIM_VERSION version ) const;
private:
YASIM_VERSION _version;
};
inline bool Version::isVersion( YASIM_VERSION version )
inline bool Version::isVersion( YASIM_VERSION version ) const
{
return _version == version;
}
inline bool Version::isVersionOrNewer( YASIM_VERSION version )
inline bool Version::isVersionOrNewer( YASIM_VERSION version ) const
{
return _version >= version;
}

42
src/FDM/YASim/Wing.cpp
View file

@ -7,47 +7,12 @@ static const float RAD2DEG = 57.2957795131;
Wing::Wing( Version * version ) :
_version(version)
{
_mirror = false;
_base[0] = _base[1] = _base[2] = 0;
_length = 0;
_chord = 0;
_taper = 0;
_sweep = 0;
_dihedral = 0;
_stall = 0;
_stallWidth = 0;
_stallPeak = 0;
_twist = 0;
_camber = 0;
_incidence = 0;
_inducedDrag = 1;
_dragScale = 1;
_liftRatio = 1;
_flap0Start = 0;
_flap0End = 0;
_flap0Lift = 0;
_flap0Drag = 0;
_flap1Start = 0;
_flap1End = 0;
_flap1Lift = 0;
_flap1Drag = 0;
_spoilerStart = 0;
_spoilerEnd = 0;
_spoilerLift = 0;
_spoilerDrag = 0;
_slatStart = 0;
_slatEnd = 0;
_slatAoA = 0;
_slatDrag = 0;
_meanChord = 0;
_wingspan = 0;
_aspectRatio = 1;
}
Wing::~Wing()
{
int i;
for(i=0; i<_surfs.size(); i++) {
for(int i=0; i<_surfs.size(); i++) {
SurfRec* s = (SurfRec*)_surfs.get(i);
delete s->surface;
delete s;
@ -57,8 +22,7 @@ Wing::~Wing()
void Wing::setIncidence(float incidence)
{
_incidence = incidence;
int i;
for(i=0; i<_surfs.size(); i++)
for(int i=0; i<_surfs.size(); i++)
((SurfRec*)_surfs.get(i))->surface->setIncidence(incidence);
}

80
src/FDM/YASim/Wing.hpp
View file

@ -106,55 +106,55 @@ private:
Vector _slatSurfs;
Vector _spoilerSurfs;
bool _mirror;
bool _mirror {false};
float _base[3];
float _length;
float _chord;
float _taper;
float _sweep;
float _dihedral;
float _base[3] {0,0,0};
float _length {0};
float _chord {0};
float _taper {1};
float _sweep {0};
float _dihedral {0};
// calculated from above
float _tip[3];
float _meanChord; // std. mean chord
float _mac; // mean aerodynamic chord length
float _macRootDistance; // y-distance of mac from root
float _macX; // x-coordinate of mac (leading edge)
float _netSpan;
float _wingspan;
float _aspectRatio;
float _tip[3] {0,0,0};
float _meanChord {0}; // std. mean chord
float _mac {0}; // mean aerodynamic chord length
float _macRootDistance {0}; // y-distance of mac from root
float _macX {0}; // x-coordinate of mac (leading edge)
float _netSpan {0};
float _wingspan {0};
float _aspectRatio {1};
float _stall;
float _stallWidth;
float _stallPeak;
float _twist;
float _camber;
float _incidence;
float _inducedDrag;
float _stall {0};
float _stallWidth {0};
float _stallPeak {0};
float _twist {0};
float _camber {0};
float _incidence {0};
float _inducedDrag {1};
float _dragScale;
float _liftRatio;
float _dragScale {1};
float _liftRatio {1};
float _flap0Start;
float _flap0End;
float _flap0Lift;
float _flap0Drag;
float _flap0Start {0};
float _flap0End {0};
float _flap0Lift {0};
float _flap0Drag {0};
float _flap1Start;
float _flap1End;
float _flap1Lift;
float _flap1Drag;
float _flap1Start {0};
float _flap1End {0};
float _flap1Lift {0};
float _flap1Drag {0};
float _spoilerStart;
float _spoilerEnd;
float _spoilerLift;
float _spoilerDrag;
float _spoilerStart {0};
float _spoilerEnd {0};
float _spoilerLift {0};
float _spoilerDrag {0};
float _slatStart;
float _slatEnd;
float _slatAoA;
float _slatDrag;
float _slatStart {0};
float _slatEnd {0};
float _slatAoA {0};
float _slatDrag {0};
Version * _version;
};

36
src/FDM/YASim/YASim.cxx
View file

@ -15,6 +15,7 @@
#include <Main/globals.hxx>
#include <Main/fg_props.hxx>
#include "yasim-common.hpp"
#include "FGFDM.hpp"
#include "Atmosphere.hpp"
#include "Math.hpp"
@ -34,22 +35,6 @@
using namespace yasim;
using std::string;
static const float YASIM_PI = 3.14159265358979323846;
static const float RAD2DEG = 180/YASIM_PI;
static const float PI2 = YASIM_PI*2;
static const float M2FT = 3.2808399;
static const float FT2M = 0.3048;
static const float MPS2KTS = 3600.0/1852.0;
static const float INHG2PA = 3386.389;
static const float SLUG2KG = 14.59390;
#if 0 // unused constant values
static const float RAD2RPM = 9.54929658551;
static const float CM2GALS = 264.172037284; // gallons/cubic meter
static const float KG2LBS = 2.20462262185;
static const float W2HP = 1.3416e-3;
#endif
YASim::YASim(double dt) :
_simTime(0)
{
@ -312,6 +297,10 @@ void YASim::copyToYASim(bool copyState)
float temp = _temp_degc->getFloatValue() + 273.15;
float dens = _density_slugft3->getFloatValue() *
SLUG2KG * M2FT*M2FT*M2FT;
Atmosphere atmo;
atmo.setDensity(dens);
atmo.setTemperature(temp);
atmo.setPressure(pressure);
// Convert and set:
Model* model = _fdm->getAirplane()->getModel();
@ -332,7 +321,7 @@ void YASim::copyToYASim(bool copyState)
Math::mmul33(s.orient, xyz2ned, s.orient);
// Velocity
float v[3];
float v[3] {0, 0, 0};
bool needCopy = false;
switch (_speed_set) {
case NED:
@ -347,15 +336,14 @@ void YASim::copyToYASim(bool copyState)
Math::tmul33(s.orient, v, v);
break;
case KNOTS:
v[0] = Atmosphere::spdFromVCAS(get_V_calibrated_kts()/MPS2KTS,
pressure, temp);
v[0] = atmo.spdFromVCAS(get_V_calibrated_kts()/MPS2KTS);
v[1] = 0;
v[2] = 0;
Math::tmul33(s.orient, v, v);
needCopy = true;
break;
case MACH:
v[0] = Atmosphere::spdFromMach(get_Mach_number(), temp);
v[0] = atmo.spdFromMach(get_Mach_number());
v[1] = 0;
v[2] = 0;
Math::tmul33(s.orient, v, v);
@ -374,21 +362,19 @@ void YASim::copyToYASim(bool copyState)
if(copyState || needCopy)
model->setState(&s);
// wind
Math::tmul33(xyz2ned, wind, wind);
model->setWind(wind);
// air
model->setAir(pressure, temp, dens);
model->setAtmosphere(atmo);
// Query a ground plane for each gear/hook/launchbar and
// write that value into the corresponding class.
_fdm->getAirplane()->getModel()->updateGround(&s);
Launchbar* l = model->getLaunchbar();
if (l)
if (l) {
l->setLaunchCmd(0.0 < _catapult_launch_cmd->getFloatValue());
}
}
// All the settables:
//

21
src/FDM/YASim/proptest.cpp
View file

@ -29,23 +29,6 @@ float fgGetFloat (const char * name, float defaultValue) { return 0; }
double fgGetDouble (const char * name, double defaultValue) { return 0; }
bool fgSetDouble (const char * name, double defaultValue) { return false; }
static const float KTS2MPS = 0.514444444444;
static const float RPM2RAD = 0.10471975512;
static const float HP2W = 745.700;
static const float FT2M = 0.3048;
static const float N2LB = 0.224809;
// static const float DEG2RAD = 0.0174532925199;
// static const float LBS2N = 4.44822;
// static const float LBS2KG = 0.45359237;
// static const float KG2LBS = 2.2046225;
// static const float CM2GALS = 264.172037284;
// static const float INHG2PA = 3386.389;
// static const float K2DEGF = 1.8;
// static const float K2DEGFOFFSET = -459.4;
// static const float CIN2CM = 1.6387064e-5;
// static const float YASIM_PI = 3.14159265358979323846;
const int COUNT = 100;
int main(int argc, char** argv)
@ -71,9 +54,7 @@ int main(int argc, char** argv)
eng->setBoost(1);
float alt = (argc > 2 ? atof(argv[2]) : 0) * FT2M;
pe->setAir(Atmosphere::getStdPressure(alt),
Atmosphere::getStdTemperature(alt),
Atmosphere::getStdDensity(alt));
pe->setStandardAtmosphere(alt);
float speed = (argc > 3 ? atof(argv[3]) : 0) * KTS2MPS;
float wind[3];

12
src/FDM/YASim/yasim-atmotest.cpp Normal file
View file

@ -0,0 +1,12 @@
#include "Atmosphere.hpp"
using namespace yasim;
int main(int argc, char** argv)
{
Atmosphere a;
if (a.test()) {
return 0;
}
return 1;
}

34
src/FDM/YASim/yasim-common.hpp Normal file
View file

@ -0,0 +1,34 @@
#ifndef _YASIM_COMMON_HPP
#define _YASIM_COMMON_HPP
namespace yasim {
static const float YASIM_PI = 3.14159265358979323846;
static const float PI2 = YASIM_PI*2;
static const float RAD2DEG = 180/YASIM_PI;
static const float DEG2RAD = YASIM_PI/180;
static const float RPM2RAD = YASIM_PI/30;
static const float KTS2MPS = 1852.0/3600.0;
static const float MPS2KTS = 3600.0/1852.0;
static const float KMH2MPS = 1/3.6;
static const float FT2M = 0.3048;
static const float M2FT = 1/FT2M;
static const float LBS2N = 4.44822;
static const float N2LB = 1/LBS2N;
static const float LBS2KG = 0.45359237;
static const float KG2LBS = 1/LBS2KG;
static const float CM2GALS = 264.172037284;
static const float HP2W = 745.700;
static const float INHG2PA = 3386.389;
static const float K2DEGF = 1.8;
static const float K2DEGFOFFSET = -459.4;
static const float CIN2CM = 1.6387064e-5;
static const float NM2FTLB = (1/(LBS2N*FT2M));
static const float SLUG2KG = 14.59390;
};
#endif // ifndef _YASIM_COMMON_HPP

24
src/FDM/YASim/yasim-test.cpp
View file

@ -7,6 +7,7 @@
#include <simgear/xml/easyxml.hxx>
#include <simgear/misc/sg_path.hxx>
#include "yasim-common.hpp"
#include "FGFDM.hpp"
#include "Atmosphere.hpp"
#include "RigidBody.hpp"
@ -26,11 +27,6 @@ float fgGetFloat (const char * name, float defaultValue) { return 0; }
double fgGetDouble (const char * name, double defaultValue = 0.0) { return 0; }
bool fgSetDouble (const char * name, double defaultValue = 0.0) { return 0; }
static const float RAD2DEG = 57.2957795131;
static const float DEG2RAD = 0.0174532925199;
/// knots 2 meters per second
static const float KTS2MPS = 0.514444444444;
enum Config
{
@ -54,7 +50,7 @@ void yasim_graph(Airplane* a, const float alt, const float kts, int cfg = CONFIG
Model* m = a->getModel();
State s;
m->setAirFromStandardAtmosphere(alt);
m->setStandardAtmosphere(alt);
switch (cfg) {
case CONFIG_APPROACH:
@ -74,7 +70,7 @@ void yasim_graph(Airplane* a, const float alt, const float kts, int cfg = CONFIG
for(int deg=-15; deg<=90; deg++) {
float aoa = deg * DEG2RAD;
Airplane::setupState(aoa, kts * KTS2MPS, 0 ,&s);
s.setupState(aoa, kts * KTS2MPS, 0);
m->getBody()->reset();
m->initIteration();
m->calcForces(&s);
@ -129,7 +125,7 @@ void yasim_drag(Airplane* a, const float aoa, const float alt, int cfg = CONFIG_
Model* m = a->getModel();
State s;
m->setAirFromStandardAtmosphere(alt);
m->setStandardAtmosphere(alt);
switch (cfg) {
case CONFIG_APPROACH:
@ -148,7 +144,7 @@ void yasim_drag(Airplane* a, const float aoa, const float alt, int cfg = CONFIG_
printf("#kts, drag\n");
for(int kts=15; kts<=150; kts++) {
Airplane::setupState(aoa, kts * KTS2MPS, 0 ,&s);
s.setupState(aoa, kts * KTS2MPS, 0);
m->getBody()->reset();
m->initIteration();
m->calcForces(&s);
@ -254,11 +250,11 @@ int main(int argc, char** argv)
float MACy = a->getWing()->getMACy();
printf(" Iterations: %d\n", a->getSolutionIterations());
printf(" Drag Coefficient: %f\n", drag);
printf(" Lift Ratio: %f\n", a->getLiftRatio());
printf(" Cruise AoA: %f deg\n", aoa);
printf(" Tail Incidence: %f deg\n", tail);
printf("Approach Elevator: %f\n\n", a->getApproachElevator());
printf(" Drag Coefficient: %.3f\n", drag);
printf(" Lift Ratio: %.3f\n", a->getLiftRatio());
printf(" Cruise AoA: %.2f deg\n", aoa);
printf(" Tail Incidence: %.2f deg\n", tail);
printf("Approach Elevator: %.3f\n\n", a->getApproachElevator());
printf(" CG: x:%.3f, y:%.3f, z:%.3f\n", cg[0], cg[1], cg[2]);
printf(" Wing MAC (*1): x:%.2f, y:%.2f, length:%.1f \n", MACx, MACy, MAC);
printf(" CG-x rel. MAC: %.3f\n", a->getCGMAC());