YASim: more output for CLI tool, e.g. wingspan, area, levers, load factor...

src/FDM/YASim/Airplane.cpp

@@ -511,7 +511,7 @@ void Airplane::compile()
     // The Wing objects
     if (_wing)
     {
-        if (baseN != 0) {
+        if (baseN != nullptr) {
             _wingsN = baseN->getChild("wing", 0, true);
             _wing->setPropertyNode(_wingsN);
         }
@@ -528,14 +528,14 @@ void Airplane::compile()
     }
     if (_tail)
     {
-        if (baseN != 0) {
+        if (baseN != nullptr) {
             _wingsN = baseN->getChild("tail", 0, true);
             _tail->setPropertyNode(_wingsN);
         }
         aeroWgt += compileWing(_tail);
     }
-    int i;
+
-    for(i=0; i<_vstabs.size(); i++)
+    for(int i=0; i<_vstabs.size(); i++)
     {
         Wing* vs = (Wing*)_vstabs.get(i);
         if (baseN != 0) {
@@ -546,32 +546,37 @@ void Airplane::compile()
     }
 
     // The fuselage(s)
-    for(i=0; i<_fuselages.size(); i++)
+    for(int i=0; i<_fuselages.size(); i++)
         aeroWgt += compileFuselage((Fuselage*)_fuselages.get(i));
 
     // Count up the absolute weight we have
     float nonAeroWgt = _ballast;
-    for(i=0; i<_thrusters.size(); i++)
+    for(int i=0; i<_thrusters.size(); i++)
         nonAeroWgt += ((ThrustRec*)_thrusters.get(i))->mass;
 
     // Rescale to the specified empty weight
     float wscale = (_emptyWeight-nonAeroWgt)/aeroWgt;
-    for(i=firstMass; i<body->numMasses(); i++) {
+    for(int i=firstMass; i<body->numMasses(); i++) {
         body->setMass(i, body->getMass(i)*wscale);
     }
-    if (_wingsN != nullptr) {
+    //if we have prop tree, give scale factor to each wing so it can export its mass to the prop tree
-        float w = _wingsN->getNode("weight", true)->getFloatValue();
+    if (baseN != nullptr) {
-        _wingsN->getNode("mass", true)->setFloatValue(w * wscale);
+        if (_wing) _wing->weight2mass(wscale);
+        if (_tail) _tail->weight2mass(wscale);
+        for(int i=0; i<_vstabs.size(); i++)
+        {
+            ((Wing*)_vstabs.get(i))->weight2mass(wscale);
+        }
     }
     // Add the thruster masses
-    for(i=0; i<_thrusters.size(); i++) {
+    for(int i=0; i<_thrusters.size(); i++) {
         ThrustRec* t = (ThrustRec*)_thrusters.get(i);
         body->addMass(t->mass, t->cg, true);
     }
 
     // Add the tanks, empty for now.
     float totalFuel = 0;
-    for(i=0; i<_tanks.size(); i++) { 
+    for(int i=0; i<_tanks.size(); i++) { 
         Tank* t = (Tank*)_tanks.get(i); 
         t->handle = body->addMass(0, t->pos);
         totalFuel += t->cap;
@@ -583,11 +588,11 @@ void Airplane::compile()
     body->recalc();
 
     // Add surfaces for the landing gear.
-    for(i=0; i<_gears.size(); i++)
+    for(int i=0; i<_gears.size(); i++)
         compileGear((GearRec*)_gears.get(i));
 
     // The Thruster objects
-    for(i=0; i<_thrusters.size(); i++) {
+    for(int i=0; i<_thrusters.size(); i++) {
         ThrustRec* tr = (ThrustRec*)_thrusters.get(i);
         tr->handle = _model.addThruster(tr->thruster);
     }
@@ -727,54 +732,52 @@ void Airplane::setupWeights(bool isApproach)
     }
 }
 
-/// load values for controls as defined in cruise/approach configuration
+void Airplane::setControlValues(const Vector& controls)
-void Airplane::loadControls(const Vector& controls)
 {
-  _controls.reset();
+    _controls.reset();
-  for(int i=0; i < controls.size(); i++) {
+    for(int i=0; i < controls.size(); i++) {
         ControlSetting* c = (ControlSetting*)controls.get(i);
-    _controls.setInput(c->propHandle, c->val);
+        _controls.setInput(c->propHandle, c->val);
-  }
+    }
-  _controls.applyControls(); 
+    _controls.applyControls(); 
 }
 
-/// Helper for solve()
 void Airplane::runConfig(Config &cfg)
 {
-  // aoa is consider to be given for approach so we calculate orientation 
+    // aoa is consider to be given for approach so we calculate orientation 
-  // only once in setApproach()
+    // only once in setApproach()
-  if (!cfg.isApproach) {
+    if (!cfg.isApproach) {
-    cfg.state.setupOrientationFromAoa(cfg.aoa);
+        cfg.state.setupOrientationFromAoa(cfg.aoa);
-  }
+    }
-  cfg.state.setupSpeedAndPosition(cfg.speed, cfg.glideAngle);
+    cfg.state.setupSpeedAndPosition(cfg.speed, cfg.glideAngle);
-  _model.setState(&cfg.state);
+    _model.setState(&cfg.state);
-  _model.setStandardAtmosphere(cfg.altitude);
+    _model.setStandardAtmosphere(cfg.altitude);
-  loadControls(cfg.controls);
+    setControlValues(cfg.controls);
   
-  // The local wind
+    // The local wind
-  float wind[3];
+    float wind[3];
-  Math::mul3(-1, cfg.state.v, wind);
+    Math::mul3(-1, cfg.state.v, wind);
-  cfg.state.globalToLocal(wind, wind);
+    cfg.state.globalToLocal(wind, wind);
-  
+    
-  setFuelFraction(cfg.fuel);
+    setFuelFraction(cfg.fuel);
-  setupWeights(cfg.isApproach);
+    setupWeights(cfg.isApproach);
-  
+    
-  // Set up the thruster parameters and iterate until the thrust
+    // Set up the thruster parameters and iterate until the thrust
-  // stabilizes.
+    // stabilizes.
-  for(int i=0; i<_thrusters.size(); i++) {
+    for(int i=0; i<_thrusters.size(); i++) {
-    Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
+        Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
-    t->setWind(wind);
+        t->setWind(wind);
-    t->setStandardAtmosphere(cfg.altitude);
+        t->setStandardAtmosphere(cfg.altitude);
-  }
+    }
-  
+    
-  stabilizeThrust();
+    stabilizeThrust();
-  updateGearState();
+    updateGearState();
-  
+    
-  // Precompute thrust in the model, and calculate aerodynamic forces
+    // Precompute thrust in the model, and calculate aerodynamic forces
-  _model.getBody()->recalc();
+    _model.getBody()->recalc();
-  _model.getBody()->reset();
+    _model.getBody()->reset();
-  _model.initIteration();
+    _model.initIteration();
-  _model.calcForces(&cfg.state);
+    _model.calcForces(&cfg.state);
 }
 /// Used only in Airplane::solve() and solveHelicopter(), not at runtime
 void Airplane::applyDragFactor(float factor)
@@ -846,6 +849,24 @@ float Airplane::normFactor(float f)
     return f;
 }
 
+///helper for Airplane::solve()
+float Airplane::_getPitch(Config &cfg)
+{
+    float tmp[3];
+    _model.getBody()->getAngularAccel(tmp);
+    cfg.state.localToGlobal(tmp, tmp);
+    return tmp[1];
+}
+
+///helper for Airplane::solve()
+float Airplane::_getLift(Config &cfg)
+{
+    float tmp[3];
+    _model.getBody()->getAccel(tmp);
+    cfg.state.localToGlobal(tmp, tmp);
+    return cfg.weight * tmp[2];
+}
+
 void Airplane::solve()
 {
     static const float ARCMIN = 0.0002909f;
@@ -859,50 +880,35 @@ void Airplane::solve()
             _failureMsg = "Solution failed to converge after 10000 iterations";
             return;
         }
-
         // Run an iteration at cruise, and extract the needed numbers:
         runConfig(_cruiseConfig);
-
         _model.getThrust(tmp);
         float thrust = tmp[0] + _cruiseConfig.weight * Math::sin(_cruiseConfig.glideAngle) * 9.81;
 
         _model.getBody()->getAccel(tmp);
         _cruiseConfig.state.localToGlobal(tmp, tmp);
         float xforce = _cruiseConfig.weight * tmp[0];
-        float clift0 = _cruiseConfig.weight * tmp[2];
+        float clift0 = _getLift(_cruiseConfig);
-
+        float pitch0 = _getPitch(_cruiseConfig);
-        _model.getBody()->getAngularAccel(tmp);
-        _cruiseConfig.state.localToGlobal(tmp, tmp);
-        float pitch0 = tmp[1];
 
         // Run an approach iteration, and do likewise
         runConfig(_approachConfig);
-
+        double apitch0 = _getPitch(_approachConfig);
-        _model.getBody()->getAngularAccel(tmp);
+        float alift = _getLift(_approachConfig);
-        _approachConfig.state.localToGlobal(tmp, tmp);
-        double apitch0 = tmp[1];
-
-        _model.getBody()->getAccel(tmp);
-        _approachConfig.state.localToGlobal(tmp, tmp);
-        float alift = _approachConfig.weight * tmp[2];
 
         // Modify the cruise AoA a bit to get a derivative
         _cruiseConfig.aoa += ARCMIN;
         runConfig(_cruiseConfig);
         _cruiseConfig.aoa -= ARCMIN;
             
-        _model.getBody()->getAccel(tmp);
+        float clift1 = _getLift(_cruiseConfig);
-        _cruiseConfig.state.localToGlobal(tmp, tmp);
-        float clift1 = _cruiseConfig.weight * tmp[2];
 
         // Do the same with the tail incidence
         _tail->setIncidence(_tailIncidence + ARCMIN);
         runConfig(_cruiseConfig);
         _tail->setIncidence(_tailIncidence);
 
-        _model.getBody()->getAngularAccel(tmp);
+        float pitch1 = _getPitch(_cruiseConfig);
-        _cruiseConfig.state.localToGlobal(tmp, tmp);
-        float pitch1 = tmp[1];
 
         // Now calculate:
         float awgt = 9.8f * _approachConfig.weight;
@@ -925,9 +931,7 @@ void Airplane::solve()
         runConfig(_approachConfig);
         _approachElevator.val -= ELEVDIDDLE;
 
-        _model.getBody()->getAngularAccel(tmp);
+        double apitch1 = _getPitch(_approachConfig);
-        _approachConfig.state.localToGlobal(tmp, tmp);
-        double apitch1 = tmp[1];
         float elevDelta = -apitch0 * (ELEVDIDDLE/(apitch1-apitch0));
 
         // Now apply the values we just computed.  Note that the
@@ -1021,4 +1025,53 @@ float Airplane::getCGMAC()
     return 0;
 }
 
+float Airplane::getWingSpan() const
+{
+    if (_wing == nullptr) return -1;
+    return  _wing->getSpan();
+}
+
+float Airplane::getWingArea() const
+{
+    if (_wing == nullptr) return -1;
+    return  _wing->getArea();
+}
+
+float Airplane::_getWingLoad(float mass) const
+{
+    if (_wing == nullptr) return -1;
+    float area =  _wing->getArea();
+    if (area == 0) return -1;
+    else return mass / area;
+}
+
+/// get x-distance between CG and 25% MAC of w
+float Airplane::_getWingLever(Wing* w) const
+{
+    if (w == nullptr) return -1;
+    float cg[3];
+    _model.getCG(cg);
+    // aerodynamic center is at 25% of MAC
+    float ac = w->getMACx() - 0.25f * w->getMACLength();
+    return ac - cg[0];
+}
+
+/// get max thrust with standard atmosphere at sea level
+float Airplane::getMaxThrust()
+{
+    float wind[3] {0,0,0};
+    float thrust[3] {0,0,0};
+    float sum[3] {0,0,0};
+    for(int i=0; i<_thrusters.size(); i++) {
+        Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
+        t->setWind(wind);
+        t->setStandardAtmosphere(0);
+        t->setThrottle(1);
+        t->stabilize();
+        t->getThrust(thrust);
+        Math::add3(thrust, sum, sum);
+    }    
+    return sum[0];
+}
+
 }; // namespace yasim

src/FDM/YASim/Airplane.hpp

@@ -37,20 +37,20 @@ public:
 
     void setPilotPos(float* pos) { Math::set3(pos, _pilotPos); }
     void getPilotPos(float* out) { Math::set3(_pilotPos, out); }
-
     void getPilotAccel(float* out);
 
-    void setEmptyWeight(float weight) {  _emptyWeight = weight; }
+    void  setEmptyWeight(float weight) {  _emptyWeight = weight; }
+    float getEmptyWeight() const {  return _emptyWeight; }
 
     Wing* getWing();
-    bool hasWing() const { return (_wing != nullptr); }
+    bool  hasWing() const { return (_wing != nullptr); }
     Wing* getTail(); 
-    void addVStab(Wing* vstab) { _vstabs.add(vstab); }
+    void  addVStab(Wing* vstab) { _vstabs.add(vstab); }
 
     void addFuselage(float* front, float* back, float width,
                      float taper=1, float mid=0.5f,
                      float cx=1, float cy=1, float cz=1, float idrag=1);
-    int addTank(float* pos, float cap, float fuelDensity);
+    int  addTank(float* pos, float cap, float fuelDensity);
     void addGear(Gear* g);
     void addHook(Hook* h) { _model.addHook(h); }
     void addLaunchbar(Launchbar* l) { _model.addLaunchbar(l); }
@@ -58,7 +58,7 @@ public:
     void addBallast(float* pos, float mass);
     void addHitch(Hitch* h) { _model.addHitch(h); }
 
-    int addWeight(float* pos, float size);
+    int  addWeight(float* pos, float size);
     void setWeight(int handle, float mass);
 
     void setApproach(float speed, float altitude, float aoa, float fuel, float gla);
@@ -83,11 +83,11 @@ public:
     
     int numTanks() const { return _tanks.size(); }
     void setFuelFraction(float frac); // 0-1, total amount of fuel
-    // get fuel in kg
+    /// get fuel in kg
     float getFuel(int tank) const { return ((Tank*)_tanks.get(tank))->fill; }
-    // set fuel in kg
+    /// set fuel in kg
     float setFuel(int tank, float fuel) { return ((Tank*)_tanks.get(tank))->fill = fuel; }
-    // get fuel density in kg/m^3
+    /// get fuel density in kg/m^3
     float getFuelDensity(int tank) const { return ((Tank*)_tanks.get(tank))->density; }
     float getTankCapacity(int tank) const { return ((Tank*)_tanks.get(tank))->cap; }
 
@@ -104,17 +104,30 @@ public:
     float getApproachElevator() const { return _approachElevator.val; }
     const char* getFailureMsg() const { return _failureMsg; }
 
-    void loadApproachControls() { loadControls(_approachConfig.controls); }
+    void setApproachControls() { setControlValues(_approachConfig.controls); }
-    void loadCruiseControls() { loadControls(_cruiseConfig.controls); }
+    void setCruiseControls() { setControlValues(_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)
     float getCGMAC(); // return c.g. x as fraction of MAC
-    // set desired range for C.G. in % of MAC, 0% = leading edge, 100% trailing edge
+    /// set desired range for C.G. in % of MAC, 0% = leading edge, 100% trailing edge
     void  setDesiredCGRangeInPercentOfMAC(float MACPercentMin, float MACPercentMax) { _cgDesiredMin = MACPercentMin; _cgDesiredMax = MACPercentMax; }
     float getCGSoftLimitXMin() const { return _cgDesiredAft; }   // get x-coordinate limit calculated from MAC and setCGRange values
     float getCGSoftLimitXMax() const { return _cgDesiredFront; } // get x-coordinate limit calculated from MAC and setCGRange values
-    void  setAutoBallast(bool allowed) { _autoBallast = allowed; }
+
+    void  setMTOW(float mtow) { _mtow = mtow; }
+    float getMTOW() const { return _mtow; }
+    float getWingSpan() const;
+    float getWingArea() const;
+    float getWingLoadEmpty() const { return _getWingLoad(_emptyWeight); };
+    float getWingLoadMTOW() const { return _getWingLoad(_mtow); };
+    /// get x-distance between CG and 25% MAC of wing
+    float getWingLever() const { return _getWingLever(_wing); };
+    /// get x-distance between CG and 25% MAC of tail
+    float getTailLever() const { return _getWingLever(_tail); };
+    float getMaxThrust();
+    float getThrust2WeightEmpty() { return getMaxThrust()/(_emptyWeight * KG2N); };
+    float getThrust2WeightMTOW() { return getMaxThrust()/(_mtow*KG2N); };
     
 private:
     struct Tank { 
@@ -161,17 +174,21 @@ private:
       float fuel {0};
       float glideAngle {0};      
       float aoa {0};
-      float altitude;
+      float altitude {0};
-      float weight;
+      float weight {0};
       State state;
       Vector controls;
     };
     Config _cruiseConfig;
     Config _approachConfig;
 
-    void loadControls(const Vector& controls);
+    /// load values for controls as defined in cruise/approach configuration
+    void setControlValues(const Vector& controls);
+    /// Helper for solve()
     void runConfig(Config &cfg);
     void solveGear();
+    float _getPitch(Config &cfg);
+    float _getLift(Config &cfg);
     void solve();
     void solveHelicopter();
     float compileWing(Wing* w);
@@ -186,11 +203,17 @@ private:
     void updateGearState();
     void setupWeights(bool isApproach);
     void calculateCGHardLimits();
+    ///calculate mass divided by area of main wing
+    float _getWingLoad(float mass) const;
+    ///calculate distance between CGx and AC of wing w
+    float _getWingLever(Wing* w) const;
     
     Model _model;
     ControlMap _controls;
 
     float _emptyWeight {0};
+    ///max take of weight
+    float _mtow {0};
     float _pilotPos[3] {0, 0, 0};
 
     Wing* _wing {nullptr};
@@ -216,14 +239,18 @@ private:
     float _tailIncidence {0};
     ControlSetting _approachElevator;
     const char* _failureMsg {0};
-    
+    /// hard limits for cg from gear position
-    float _cgMax {-1e6};         // hard limits for cg from gear position
+    float _cgMax {-1e6};         
-    float _cgMin {1e6};          // hard limits for cg from gear position
+    /// hard limits for cg from gear position
-    float _cgDesiredMax {0.3f};  // desired cg max in %MAC from config
+    float _cgMin {1e6};          
-    float _cgDesiredMin {0.25f}; // desired cg min in %MAC from config
+    /// desired cg max in %MAC from config
-    float _cgDesiredFront {0};   // calculated desired cg x max
+    float _cgDesiredMax {0.3f};  
-    float _cgDesiredAft {0};     // calculated desired cg x min 
+    /// desired cg min in %MAC from config
-    bool _autoBallast = false;
+    float _cgDesiredMin {0.25f}; 
+    /// calculated desired cg x max
+    float _cgDesiredFront {0};   
+    /// calculated desired cg x min 
+    float _cgDesiredAft {0};     
 };
 
 }; // namespace yasim

src/FDM/YASim/FGFDM.cpp

@@ -285,7 +285,10 @@ void FGFDM::parseAirplane(const XMLAttributes* a)
         SG_LOG(SG_FLIGHT, SG_DEV_ALERT, "This aircraft does not use the latest yasim configuration version.");
     }
     _airplane.setDesiredCGRangeInPercentOfMAC(attrf(a, "cg-min", 0.25f), attrf(a, "cg-max", 0.3f)); 
-    if (attrb(a, "auto-ballast")) { _airplane.setAutoBallast(true); }
+    
+    if (a->hasAttribute("mtow-lbs")) { _airplane.setMTOW(attrf(a, "mtow-lbs") * LBS2KG); }
+    else if (a->hasAttribute("mtow-kg")) { _airplane.setMTOW(attrf(a, "mtow-kg")); }
+    
 }
 
 void FGFDM::parseApproachCruise(const XMLAttributes* a, const char* name)
@@ -566,6 +569,8 @@ void FGFDM::parseWing(const XMLAttributes* a, const char* type, Airplane* airpla
     } 
     else if (!strcmp(type, "hstab")) {
         w = airplane->getTail();
+        if (a->hasAttribute("incidence-max")) w->setIncidenceMax(attrf(a, "incidence-max") * DEG2RAD);
+        if (a->hasAttribute("incidence-min")) w->setIncidenceMin(attrf(a, "incidence-min") * DEG2RAD);
     } else {
         w = new Wing(airplane, mirror);
     }
@@ -816,22 +821,22 @@ void FGFDM::parsePropeller(const XMLAttributes* a)
     prop->setStops(fine_stop, coarse_stop);
 
     if(a->hasAttribute("takeoff-power")) {
-	float power0 = attrf(a, "takeoff-power") * HP2W;
+        float power0 = attrf(a, "takeoff-power") * HP2W;
-	float omega0 = attrf(a, "takeoff-rpm") * RPM2RAD;
+        float omega0 = attrf(a, "takeoff-rpm") * RPM2RAD;
-	prop->setTakeoff(omega0, power0);
+        prop->setTakeoff(omega0, power0);
     }
 
     if(a->hasAttribute("max-rpm")) {
-	float max = attrf(a, "max-rpm") * RPM2RAD;
+        float max = attrf(a, "max-rpm") * RPM2RAD;
-	float min = attrf(a, "min-rpm") * RPM2RAD;
+        float min = attrf(a, "min-rpm") * RPM2RAD;
-	thruster->setVariableProp(min, max);
+        thruster->setVariableProp(min, max);
     }
 
     if(attrb(a, "contra"))
         thruster->setContraPair(true);
 
     if(a->hasAttribute("manual-pitch")) {
-	prop->setManualPitch();
+        prop->setManualPitch();
     }
 
     thruster->setGearRatio(attrf(a, "gear-ratio", 1));

src/FDM/YASim/Model.hpp

@@ -29,6 +29,7 @@ public:
     virtual ~Model();
 
     RigidBody* getBody() { return &_body; }
+    void getCG(float* cg) const { return _body.getCG(cg); }
     Integrator* getIntegrator() { return &_integrator; }
 
     void setTurbulence(Turbulence* turb) { _turb = turb; }

src/FDM/YASim/RigidBody.hpp

@@ -59,18 +59,19 @@ public:
     // regenerate its internal tables.  This step is expensive, so
     // it's exposed to the client who can amortize the call across
     // multiple changes. see also _recalcStatic() 
+    /// calculate the total mass, centre of gravity and inertia tensor
     void recalc();
 
-    // Resets the current force/torque parameters to zero.
+    /// Resets the current force/torque parameters to zero.
     void reset();
 
-    // Applies a force at the center of gravity.
+    /// Applies a force at the center of gravity.
     void addForce(const float* force) { Math::add3(_force, force, _force); }
     
-    // Applies a force at the specified position.
+    /// Applies a force at the specified position.
     void addForce(const float* pos, const float* force); 
     
-    // Adds a torque with the specified axis and magnitude
+    /// Adds a torque with the specified axis and magnitude
     void addTorque(const float* torque) { Math::add3(_torque, torque, _torque); }
 
     // Sets the rotation rate of the body (about its c.g.) within the

src/FDM/YASim/Wing.cpp

@@ -47,24 +47,30 @@ int Wing::addWingSection(float* base, float chord, float wingLength, float taper
     ws->_twist = twist;
     ws->_camber = camber;
     ws->_inducedDrag = idrag;
+    ws->_id = _sections.add(ws);
     ws->calculateGeometry();
-    int idx = _sections.add(ws);    
     // first / only section
-    if (idx == 0) {
+    if (ws->_id == 0) {
         _mac = ws->getMAC();
-        _wingspan = ws->_wingspan;
+        _netSpan = ws->_sectionSpan;
         _area = ws->getArea();
         _meanChord = ws->_meanChord;
+        _sweepLEMin = _sweepLEMax = ws->calculateSweepAngleLeadingEdge();
     }
     // append section: Calculate wing MAC from MACs of section and prev wing
     else {
         _mac = Wing::calculateMAC(_mac, ws->getMAC());
-        _wingspan += ws->_wingspan;
+        _netSpan += ws->_sectionSpan;
         _area += ws->getArea();
         _meanChord = _meanChord * ws->_meanChord * 0.5f;
+        float s = ws->calculateSweepAngleLeadingEdge();
+        if (_sweepLEMax < s) _sweepLEMax = s;
+        if (_sweepLEMin > s) _sweepLEMin = s;
     }
     _chord2float(ws->_tipChord, _tip);    
-    return idx;
+    _wingSpan = 2 * _tip[1];
+    _taper = ws->_tipChord.length / ((WingSection*)_sections.get(0))->_rootChord.length;
+    return ws->_id;
 }
 
 Chord Wing::_float2chord(float* pos, float lenght)
@@ -135,19 +141,22 @@ void Wing::WingSection::calculateTipChord() {
 void Wing::WingSection::calculateSpan()
 {
     // wingspan in y-direction (not for vstab)
-    _wingspan = Math::abs(2*_tipChord.y);
+    _sectionSpan = Math::abs(_rootChord.y - _tipChord.y);
-    _aspectRatio = _wingspan / _meanChord;  
 }
 
 void Wing::WingSection::calculateMAC()
-{
+{    
-    //FIXME call static method, use absolute y values
-    
     // http://www.nasascale.org/p2/wp-content/uploads/mac-calculator.htm
-    const float commonFactor = _rootChord.length*(0.5+_taper)/(3*_rootChord.length*(1+_taper));
+    //const float commonFactor = _rootChord.length*(0.5+_taper)/(3*_rootChord.length*(1+_taper));
+    
+    const float commonFactor = (0.5+_taper)/(3*(1+_taper));
     _mac.length = _rootChord.length-(2*_rootChord.length*(1-_taper)*commonFactor);
+    // y distance to root chord
     _mac.y = Math::abs(2*(_tipChord.y-_rootChord.y))*commonFactor;
+    // MAC leading edge x = midpoint + half MAC length
     _mac.x = _rootChord.x-Math::tan(_sweepAngleCenterLine)*_mac.y + _mac.length/2;
+    _mac.z = _rootChord.z+Math::tan(_dihedral)*_mac.y;
+    //add root y to get aircraft coordinates
     _mac.y += _rootChord.y;
 }
 
@@ -169,16 +178,20 @@ void Wing::WingSection::setIncidence(float incidence)
         ((SurfRec*)_surfs.get(i))->surface->setIncidence(incidence + _sectionIncidence);
 }
 
-// root and tip (x,y) coordinates are on leading edge
+/// root and tip (x,y) coordinates are on leading edge
-Chord Wing::calculateMAC(Chord root, Chord tip)
+Chord Wing::calculateMAC(const Chord root, const Chord tip)
 {
     assert(root.length > 0);
+    // http://www.nasascale.org/p2/wp-content/uploads/mac-calculator.htm
     //taper = tip.length / root.length;
     const float commonFactor = (root.length*0.5+tip.length)/(3*(root.length+tip.length));
     Chord m;
     m.length = root.length-(2*(root.length - tip.length)*commonFactor);
-    m.y = Math::abs(2*(tip.y - root.y))*commonFactor + root.y; 
+    float dy = tip.y - root.y;
-    m.x = root.x - (root.x - tip.x)*(root.y - m.y)/(root.y - tip.y);
+    m.y = Math::abs(2*dy)*commonFactor; 
+    m.z = root.z + m.y*(tip.z-root.z)/dy;
+    m.x = root.x - m.y*(root.x - tip.x)/dy;
+    m.y += root.y;
     return m;
 }
 
@@ -344,6 +357,12 @@ void Wing::compile()
     writeInfoToProptree();
 }
 
+float Wing::getArea() const 
+{ 
+    if (_mirror) return 2 * _area; 
+    else return _area;    
+};
+
 void Wing::multiplyLiftRatio(float factor)
 {
     WingSection* ws;
@@ -366,6 +385,11 @@ void Wing::multiplyDragCoefficient(float factor)
 
 void Wing::setIncidence(float incidence)
 {
+    if (incidence < _incidenceMin || incidence > _incidenceMax)
+    {
+        fprintf(stderr, "YASim: cannot set incidence, parameter out of range.");
+        return;
+    }
     WingSection* ws;
     for (int section=0; section < _sections.size(); section++) 
     {        
@@ -472,24 +496,22 @@ void Wing::writeInfoToProptree()
     if (_wingN == nullptr) 
         return;
     WingSection* ws = (WingSection*)_sections.get(0);
-    float chord = ws->_rootChord.length;
+    _wingN->getNode("root-chord", true)->setFloatValue(ws->_rootChord.length);
-    ws = (WingSection*)_sections.get(_sections.size()-1);
-    float taper = ws->_rootChord.length * ws->_taper;
     _wingN->getNode("tip-x", true)->setFloatValue(_tip[0]);
     _wingN->getNode("tip-y", true)->setFloatValue(_tip[1]);
     _wingN->getNode("tip-z", true)->setFloatValue(_tip[2]);
     _wingN->getNode("base-x", true)->setFloatValue(_base[0]);
     _wingN->getNode("base-y", true)->setFloatValue(_base[1]);
     _wingN->getNode("base-z", true)->setFloatValue(_base[2]);
-    _wingN->getNode("chord", true)->setFloatValue(chord);
+    _wingN->getNode("taper", true)->setFloatValue(getTaper());
-    _wingN->getNode("taper", true)->setFloatValue(taper);
+    _wingN->getNode("wing-span", true)->setFloatValue(getSpan());
-    _wingN->getNode("wing-span", true)->setFloatValue(_wingspan);
+    _wingN->getNode("wing-area", true)->setFloatValue(getArea());
-    _wingN->getNode("wing-area", true)->setFloatValue(_wingspan*_meanChord);
+    _wingN->getNode("aspect-ratio", true)->setFloatValue(getAspectRatio());
-    _wingN->getNode("aspect-ratio", true)->setFloatValue(_aspectRatio);
+    _wingN->getNode("standard-mean-chord", true)->setFloatValue(getSMC());
-    _wingN->getNode("standard-mean-chord", true)->setFloatValue(_meanChord);
+    _wingN->getNode("mac", true)->setFloatValue(getMACLength());
-    _wingN->getNode("mac", true)->setFloatValue(_mac.length);
+    _wingN->getNode("mac-x", true)->setFloatValue(getMACx());
-    _wingN->getNode("mac-x", true)->setFloatValue(_mac.x);
+    _wingN->getNode("mac-y", true)->setFloatValue(getMACy());
-    _wingN->getNode("mac-y", true)->setFloatValue(_mac.y);
+    _wingN->getNode("mac-z", true)->setFloatValue(getMACz());
 
     float wgt = 0;
     float dragSum = 0;
@@ -510,6 +532,14 @@ void Wing::writeInfoToProptree()
     _wingN->getNode("drag", true)->setFloatValue(dragSum);
 }
 
+void Wing::weight2mass(float scale)
+{
+    if (_wingN == nullptr) 
+        return;
+    float wgt = _wingN->getNode("weight", true)->getFloatValue();
+    _wingN->getNode("mass", true)->setFloatValue(wgt * scale);
+}
+
 // estimate a mass distibution and add masses to the model
 // they will be scaled to match total mass of aircraft later
 float Wing::updateModel(Model* model) 

src/FDM/YASim/Wing.hpp

@@ -1,6 +1,7 @@
 #ifndef _WING_HPP
 #define _WING_HPP
 
+#include "yasim-common.hpp"
 #include "Vector.hpp"
 #include "Version.hpp"
 #include "Math.hpp"
@@ -50,11 +51,12 @@ class Wing {
     
     struct WingSection {
         int _id;
+        /// length and midpoint (not leading edge!) of root chord
         Chord _rootChord;
-        // length is distance from base to tip, not wing span
+        /// length is distance from base to tip, not wing span
         float _length {0};
         float _taper {1};
-        // sweep of center line, not leading edge!
+        /// sweep of center line, not leading edge!
         float _sweepAngleCenterLine {0}; 
         float _dihedral {0};
         float _twist {0};
@@ -63,7 +65,7 @@ class Wing {
 
         StallParams _stallParams;
 
-        //fixed incidence of section as given in config XML
+        ///fixed incidence of section as given in config XML
         float _sectionIncidence {0};
         float _dragScale {1};
         float _liftRatio {1};
@@ -74,13 +76,14 @@ class Wing {
         float _orient[9];
         float _rightOrient[9];
         Chord _tipChord;
-        float _meanChord {0}; // std. mean chord
+        /// std. mean chord
-        Chord _mac; // mean aerodynamic chord (x,y) leading edge
+        float _meanChord {0}; 
-        float _wingspan {0};
+        /// mean aerodynamic chord, (x,y) leading edge!
-        float _aspectRatio {1};
+        Chord _mac; 
-        // all SurfRec of this wing
+        float _sectionSpan {0};
+        /// all SurfRec of this wing
         Vector _surfs;      
-        // surfaces having a certain type of flap (flap, slat, spoiler)
+        /// surfaces having a certain type of flap (flap, slat, spoiler)
         Vector _flapSurfs[sizeof(WingFlaps)];
         
         void calculateGeometry();
@@ -96,7 +99,7 @@ class Wing {
 
         // valid only after Wing::compile() was called
         Chord getMAC() const { return _mac; };
-        float getArea() const { return _wingspan*_meanChord; };
+        float getArea() const { return _sectionSpan*_meanChord; };
         void setDragCoefficient(float scale);
         void multiplyDragCoefficient(float factor);
         // The ratio of force along the Z (lift) direction of each wing
@@ -116,19 +119,29 @@ class Wing {
     Version* _version;
     bool _mirror {false};
     Vector _sections;
-    Chord _mac;
+    // midpoint (not leading edge!) of wing root chord
     float _base[3] {0,0,0};
+    // midpoint (not leading edge!) of wing tip
     float _tip[3] {0,0,0};
-    float _wingspan {0};
+    float _netSpan {0};
+    float _wingSpan {0};
     float _area {0};
     float _aspectRatio {0};
     float _meanChord {0};
+    Chord _mac;
+    float _taper {1};
     float _incidence {0};
+    float _incidenceMin {-20*DEG2RAD};
+    float _incidenceMax {20*DEG2RAD};
     float _weight {0};
+    float _sweepLEMin {0};
+    float _sweepLEMax {0};
 
     //-- private methods
-    Chord _float2chord(float* pos, float lenght = 0);
+    //copy float[3] to chord x,y,z
-    void _chord2float(Chord c, float* pos);
+    static Chord _float2chord(float* pos, float lenght = 0);
+    //copy chord x,y,z to float[3]
+    static void _chord2float(Chord c, float* pos);
     void interp(const float* v1, const float* v2, const float frac, float* out);
     void writeInfoToProptree();
     
@@ -139,7 +152,7 @@ public:
     int addWingSection(float* base, float chord, float wingLength, 
         float taper = 1, float sweep = 0, float dihedral = 0, float twist = 0, float camber = 0, float idrag = 1, float incidence = 0);
 
-    static Chord calculateMAC(Chord root, Chord tip);    
+    static Chord calculateMAC(const Chord root, const Chord tip);    
 
     void setFlapParams(int section, WingFlaps type, FlapParams fp);
     void setSectionDrag(int section, float pdrag);
@@ -150,17 +163,24 @@ public:
     void multiplyLiftRatio(float factor);
     void multiplyDragCoefficient(float factor);
     void setIncidence(float incidence);
-
+    void setIncidenceMin(float min) { _incidenceMin = min; };
+    void setIncidenceMax(float max) { _incidenceMax = max; };
+    void weight2mass(float scale);
+    
     bool isMirrored() const { return _mirror; };
     void getBase(float* base) const { Math::set3(_base, base); };
     void getTip(float* tip) const { Math::set3(_tip, tip); };
-    float getSpan() const { return _wingspan; };
+    float getSpan() const { return _wingSpan; };
-    float getArea() const { return _area; };
+    float getTaper() const { return _taper; };
+    float getArea() const;
     float getAspectRatio() const { return _aspectRatio; };
     float getSMC() const { return _meanChord; };
     float getMACLength() const { return _mac.length; }; // get length of MAC
     float getMACx() const { return _mac.x; }; // get x-coord of MAC leading edge 
     float getMACy() const { return _mac.y; }; // get y-coord of MAC leading edge 
+    float getMACz() const { return _mac.z; }; // get y-coord of MAC leading edge 
+    float getSweepLEMin() const { return _sweepLEMin; }; //min sweep angle of leading edge
+    float getSweepLEMax() const { return _sweepLEMax; }; //max sweep angle of leading edge
    
 //-----------------------------    
     // propergate the control axes value for the sub-surfaces

src/FDM/YASim/yasim-common.hpp

@@ -20,6 +20,7 @@ namespace yasim {
 
     static const float LBS2N = 4.44822f;
     static const float N2LB = 1/LBS2N;
+    static const float KG2N = 9.81f;
     static const float LBS2KG = 0.45359237f;
     static const float KG2LBS = 1/LBS2KG;
     static const float CM2GALS = 264.172037284f;

src/FDM/YASim/yasim-test.cpp

@@ -54,10 +54,10 @@ void yasim_graph(Airplane* a, const float alt, const float kts, int cfg = CONFIG
 
   switch (cfg) {
     case CONFIG_APPROACH:
-      a->loadApproachControls();
+      a->setApproachControls();
       break;
     case CONFIG_CRUISE:
-      a->loadCruiseControls();
+      a->setCruiseControls();
       break;
     case CONFIG_NONE:
       break;
@@ -129,10 +129,10 @@ void yasim_drag(Airplane* a, const float aoa, const float alt, int cfg = CONFIG_
   
   switch (cfg) {
     case CONFIG_APPROACH:
-      a->loadApproachControls();
+      a->setApproachControls();
       break;
     case CONFIG_CRUISE:
-      a->loadCruiseControls();
+      a->setCruiseControls();
       break;
     case CONFIG_NONE:
       break;
@@ -237,7 +237,7 @@ int main(int argc, char** argv)
     printf("= YASim solution results =\n");
     printf("==========================\n");
     float aoa = a->getCruiseAoA() * RAD2DEG;
-    float tail = -1 * a->getTailIncidence() * RAD2DEG;
+    float tailIncidence = -1 * a->getTailIncidence() * RAD2DEG;
     float drag = 1000 * a->getDragCoefficient();
     float cg[3];
     a->getModel()->getBody()->getCG(cg);
@@ -246,28 +246,49 @@ int main(int argc, char** argv)
     float SI_inertia[9];
     a->getModel()->getBody()->getInertiaMatrix(SI_inertia);
     float MAC = 0, MACx = 0, MACy = 0;
+    float sweepMin = 0, sweepMax = 0;
     Wing* wing {nullptr};
+    Wing* tail {nullptr};
     if (a->hasWing()) {
       wing = a->getWing();
-      MAC = wing->getMACLength();
+      tail = a->getTail();
-      MACx = wing->getMACx();
-      MACy = wing->getMACy();
     }
     printf("       Iterations: %d\n", a->getSolutionIterations());
     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("   Tail Incidence: %.2f deg\n", tailIncidence);
     printf("Approach Elevator: %.3f\n\n", a->getApproachElevator());
     printf("               CG: x:%.3f, y:%.3f, z:%.3f\n", cg[0], cg[1], cg[2]);
     if (wing) {
+      MAC = wing->getMACLength();
+      MACx = wing->getMACx();
+      MACy = wing->getMACy();
+      sweepMin = wing->getSweepLEMin() * RAD2DEG;
+      sweepMax = wing->getSweepLEMax() * RAD2DEG;
       printf("         Wing MAC: (x:%.2f, y:%.2f), length:%.1f \n", MACx, MACy, MAC);
-      printf("  hard limit CG-x: %.3f \n", a->getCGHardLimitXMax());
+      printf("  hard limit CG-x: %.3f m\n", a->getCGHardLimitXMax());
-      printf("  soft limit CG-x: %.3f \n", a->getCGSoftLimitXMax());
+      printf("  soft limit CG-x: %.3f m\n", a->getCGSoftLimitXMax());
-      printf("             CG-x: %.3f \n", cg[0]);
+      printf("             CG-x: %.3f m\n", cg[0]);
       printf("    CG-x rel. MAC: %3.0f%%\n", a->getCGMAC()*100);
-      printf("  soft limit CG-x: %.3f \n", a->getCGSoftLimitXMin());
+      printf("  soft limit CG-x: %.3f m\n", a->getCGSoftLimitXMin());
-      printf("  hard limit CG-x: %.3f \n", a->getCGHardLimitXMin());
+      printf("  hard limit CG-x: %.3f m\n", a->getCGHardLimitXMin());
+      printf("\n");
+      printf("        wing span: %.2f m\n", a->getWingSpan());
+      printf(" sweep lead. edge: %.1f .. %.1f deg\n", sweepMin, sweepMax);
+      printf("        wing area: %.2f m²\n", a->getWingArea());
+      printf("  wing load empty: %.2f kg/m² (Empty %.0f kg)\n", a->getWingLoadEmpty(), a->getEmptyWeight());
+      printf("   wing load MTOW: %.2f kg/m² (MTOW  %.0f kg)\n", a->getWingLoadMTOW(), a->getMTOW());
+      printf("\n");
+      printf("        tail span: %.3f m\n", tail->getSpan());
+      printf("        tail area: %.3f m²\n", tail->getArea());
+      printf("\n");
+      printf("       wing lever: %.3f m\n", a->getWingLever());
+      printf("       tail lever: %.3f m\n", a->getTailLever());
+      printf("\n");
+      printf("       max thrust: %.2f kN\n", a->getMaxThrust()/1000);
+      printf("     thrust/empty: %.2f\n", a->getThrust2WeightEmpty());
+      printf("      thrust/mtow: %.2f\n", a->getThrust2WeightMTOW());
     }
     printf("\nInertia tensor [kg*m^2], origo at CG:\n\n");
     printf("  %7.0f, %7.0f, %7.0f\n", SI_inertia[0], SI_inertia[1], SI_inertia[2]);