YASim: add MAC (mean aerodynamic chord) and c.g. calculations.

src/FDM/YASim/Airplane.cpp

@@ -57,8 +57,10 @@ Airplane::Airplane()
 
     _failureMsg = 0;
     _wingsN = 0;
-    _cgMaxX = -1e6;
+    _cgMax = -1e6;
-    _cgMinX = 1e6;
+    _cgMin = 1e6;
+    _cgDesiredMax = 0.33f; // FIXME find reasonable default value
+    _cgDesiredMin = 0.1f; // FIXME find reasonable default value
 }
 
 Airplane::~Airplane()
@@ -115,7 +117,7 @@ void Airplane::calcFuelWeights()
     }
 }
 
-void Airplane::getPilotAccel(float* out)
+const void Airplane::getPilotAccel(float* out)
 {
     State* s = _model.getState();
 
@@ -125,14 +127,13 @@ void Airplane::getPilotAccel(float* out)
     Math::vmul33(s->orient, out, out);
     out[0] = -out[0];
 
-    // The regular acceleration
+    float acceleration[3];
-    float tmp[3];
     // Convert to aircraft coordinates
-    Math::vmul33(s->orient, s->acc, tmp);
+    Math::vmul33(s->orient, s->acc, acceleration);
-    tmp[1] = -tmp[1];
+    acceleration[1] = -acceleration[1];
-    tmp[2] = -tmp[2];
+    acceleration[2] = -acceleration[2];
 
-    Math::add3(tmp, out, out);
+    Math::add3(acceleration, out, out);
 
     // FIXME: rotational & centripetal acceleration needed
 }
@@ -240,10 +241,6 @@ void Airplane::addGear(Gear* gear)
     g->gear = gear;
     g->surf = 0;
     _gears.add(g);
-    float pos[3];
-    g->gear->getPosition(pos);
-    if (pos[0] > _cgMaxX) _cgMaxX = pos[0];
-    if (pos[0] < _cgMinX) _cgMinX = pos[0];
 }
 
 void Airplane::addThruster(Thruster* thruster, float mass, float* cg)
@@ -371,6 +368,9 @@ float Airplane::compileWing(Wing* w)
       _wingsN->getNode("wing-area", true)->setFloatValue(w->getArea());
       _wingsN->getNode("aspect-ratio", true)->setFloatValue(w->getAspectRatio());
       _wingsN->getNode("standard-mean-chord", true)->setFloatValue(w->getSMC());
+      _wingsN->getNode("mac", true)->setFloatValue(w->getMAC());
+      _wingsN->getNode("mac-x", true)->setFloatValue(w->getMACx());
+      _wingsN->getNode("mac-y", true)->setFloatValue(w->getMACy());
     }
 
     float wgt = 0;
@@ -582,7 +582,6 @@ void Airplane::compile()
     RigidBody* body = _model.getBody();
     int firstMass = body->numMasses();
     SGPropertyNode_ptr baseN = fgGetNode("/fdm/yasim/model/wings", true);
-    SGPropertyNode_ptr n;
 
     // Generate the point masses for the plane.  Just use unitless
     // numbers for a first pass, then go back through and rescale to
@@ -594,6 +593,15 @@ void Airplane::compile()
     {
       if (baseN != 0) _wingsN = baseN->getChild("wing", 0, true);
       aeroWgt += compileWing(_wing);
+      
+      // convert % to absolute x coordinates
+      _cgDesiredFront = _wing->getMACx() - _wing->getMAC()*_cgDesiredMin;
+      _cgDesiredAft = _wing->getMACx() - _wing->getMAC()*_cgDesiredMax;
+      if (baseN != 0) {
+        SGPropertyNode_ptr n = fgGetNode("/fdm/yasim/model", true);
+        n->getNode("cg-range-front", true)->setFloatValue(_cgDesiredFront);
+        n->getNode("cg-range-aft", true)->setFloatValue(_cgDesiredAft);
+      }
     }
     if (_tail)
     {
@@ -649,25 +657,25 @@ void Airplane::compile()
         tr->handle = _model.addThruster(tr->thruster);
     }
     
-    // Ground effect
-    // If a double tapered wing is modelled with wing and mstab, wing must 
-    // be outboard to get correct wingspan.
     if(_wing) {
-        float gepos[3];
+        // Ground effect
+        // If a double tapered wing is modelled with wing and mstab, wing must 
+        // be outboard to get correct wingspan.
+        float pos[3];
         float gespan = 0;
         gespan = _wing->getSpan();
-        _wing->getBase(gepos);
+        _wing->getBase(pos);
         if(!isVersionOrNewer( Version::YASIM_VERSION_2017_2 )) {
           //old code
           //float span = _length * Math::cos(_sweep) * Math::cos(_dihedral);
           //span = 2*(span + Math::abs(_base[2]));
-          gespan -= 2*gepos[1]; // cut away base (y-distance)
+          gespan -= 2*pos[1]; // cut away base (y-distance)
-          gespan += 2*Math::abs(gepos[2]); // add (wrong) z-distance
+          gespan += 2*Math::abs(pos[2]); // add (wrong) z-distance
         }
         if (baseN != 0)
           baseN->getChild("wing", 0)->getNode("gnd-eff-span", true)->setFloatValue(gespan);
         // where does the hard coded factor 0.15 come from?
-        _model.setGroundEffect(gepos, gespan, 0.15f);
+        _model.setGroundEffect(pos, gespan, 0.15f);
     }
     
     // solve function below resets failure message
@@ -675,6 +683,8 @@ void Airplane::compile()
     if (_failureMsg) return;
 
     solveGear();
+    calculateCGHardLimits();
+    
     if(_wing && _tail) solve();
     else
     {
@@ -742,11 +752,24 @@ void Airplane::solveGear()
     }
 }
 
+void Airplane::calculateCGHardLimits()
+{
+    _cgMax = -1e6;
+    _cgMin = 1e6;
+    for (int i = 0; i < _gears.size(); i++) {
+        GearRec* gr = (GearRec*)_gears.get(i);
+        float pos[3];
+        gr->gear->getPosition(pos);
+        if (pos[0] > _cgMax) _cgMax = pos[0];
+        if (pos[0] < _cgMin) _cgMin = pos[0];
+    }
+}
+
 void Airplane::initEngines()
 {
     for(int i=0; i<_thrusters.size(); i++) {
         ThrustRec* tr = (ThrustRec*)_thrusters.get(i);
-	tr->thruster->init();
+        tr->thruster->init();
     }
 }
 
@@ -835,7 +858,6 @@ void Airplane::runApproach()
     Math::vmul33(_approachState.orient, wind, wind);
     
     setFuelFraction(_approachFuel);
-
     setupWeights(true);
 
     // Run the thrusters until they get to a stable setting.  FIXME:
@@ -1094,4 +1116,11 @@ void Airplane::solveHelicopter()
     
 }
 
+const float Airplane::getCGMAC()
+{ 
+    float cg[3];
+    _model.getBody()->getCG(cg);
+    return (_wing->getMACx() - cg[0]) / _wing->getMAC();
+}
+
 }; // namespace yasim

src/FDM/YASim/Airplane.hpp

@@ -33,11 +33,12 @@ public:
     void setPilotPos(float* pos) { Math::set3(pos, _pilotPos); }
     void getPilotPos(float* out) { Math::set3(_pilotPos, out); }
 
-    void getPilotAccel(float* out);
+    const void getPilotAccel(float* out);
 
     void setEmptyWeight(float weight) {  _emptyWeight = weight; }
 
     void setWing(Wing* wing) { _wing = wing; }
+    Wing* getWing() { return _wing; }
     void setTail(Wing* tail) { _tail = tail; }
     void addVStab(Wing* vstab) { _vstabs.add(vstab); }
 
@@ -64,56 +65,62 @@ public:
 
     void addSolutionWeight(bool approach, int idx, float wgt);
 
-    int numGear() { return _gears.size(); }
+    const int numGear() { return _gears.size(); }
     Gear* getGear(int g) { return ((GearRec*)_gears.get(g))->gear; }
     Hook* getHook() { return _model.getHook(); }
-    int numHitches() { return _hitches.size(); }
+    const int numHitches() { return _hitches.size(); }
     Hitch* getHitch(int h);
     Rotorgear* getRotorgear() { return _model.getRotorgear(); }
     Launchbar* getLaunchbar() { return _model.getLaunchbar(); }
 
-    int numThrusters() { return _thrusters.size(); }
+    const int numThrusters() { return _thrusters.size(); }
     Thruster* getThruster(int n) {
         return ((ThrustRec*)_thrusters.get(n))->thruster; }
     
-    int numTanks() { return _tanks.size(); }
+    const int numTanks() { return _tanks.size(); }
     void setFuelFraction(float frac); // 0-1, total amount of fuel
     // get fuel in kg
-    float getFuel(int tank) { return ((Tank*)_tanks.get(tank))->fill; }
+    const float getFuel(int tank) { return ((Tank*)_tanks.get(tank))->fill; }
     // set fuel in kg
     float setFuel(int tank, float fuel) { return ((Tank*)_tanks.get(tank))->fill = fuel; }
     // get fuel density in kg/m^3
-    float getFuelDensity(int tank) { return ((Tank*)_tanks.get(tank))->density; }
+    const float getFuelDensity(int tank) { return ((Tank*)_tanks.get(tank))->density; }
-    float getTankCapacity(int tank) { return ((Tank*)_tanks.get(tank))->cap; }
+    const float getTankCapacity(int tank) { return ((Tank*)_tanks.get(tank))->cap; }
 
     void compile(); // generate point masses & such, then solve
     void initEngines();
     void stabilizeThrust();
 
     // Solution output values
-    int getSolutionIterations() { return _solutionIterations; }
+    const int getSolutionIterations() { return _solutionIterations; }
-    float getDragCoefficient() { return _dragFactor; }
+    const float getDragCoefficient() { return _dragFactor; }
-    float getLiftRatio() { return _liftRatio; }
+    const float getLiftRatio() { return _liftRatio; }
-    float getCruiseAoA() { return _cruiseAoA; }
+    const float getCruiseAoA() { return _cruiseAoA; }
-    float getTailIncidence() { return _tailIncidence; }
+    const float getTailIncidence() { return _tailIncidence; }
-    float getApproachElevator() { return _approachElevator.val; }
+    const float getApproachElevator() { return _approachElevator.val; }
     const char* getFailureMsg() { 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); }
     
-    float getCGMinX() { return _cgMinX; }
+    const float getCGHardLimitXMin() { return _cgMin; } // get min x-coordinate for c.g (from main gear)
-    float getCGMaxX() { return _cgMaxX; }
+    const float getCGHardLimitXMax() { return _cgMax; } // get max x-coordinate for c.g (from nose gear)
+    const 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
+    void  setDesiredCGRangeInPercentOfMAC(float MACPercentMin, float MACPercentMax) { _cgDesiredMin = MACPercentMin; _cgDesiredMax = MACPercentMax; }
+    const float getCGSoftLimitXMin() { return _cgDesiredAft; }   // get x-coordinate limit calculated from MAC and setCGRange values
+    const float getCGSoftLimitXMax() { return _cgDesiredFront; } // get x-coordinate limit calculated from MAC and setCGRange values
+    void  setAutoBallast(bool allowed) { _autoBallast = allowed; }
     
 private:
     struct Tank { float pos[3]; float cap; float fill;
-	          float density; int handle; };
+        float density; int handle; };
     struct Fuselage { float front[3], back[3], width, taper, mid, _cx, _cy, _cz, _idrag;
-		      Vector surfs; };
+        Vector surfs; };
     struct GearRec { Gear* gear; Surface* surf; float wgt; };
     struct ThrustRec { Thruster* thruster;
-	               int handle; float cg[3]; float mass; };
+        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; };
@@ -136,6 +143,7 @@ private:
     float normFactor(float f);
     void updateGearState();
     void setupWeights(bool isApproach);
+    void calculateCGHardLimits();
     
     Model _model;
     ControlMap _controls;
@@ -187,9 +195,13 @@ private:
     Control _approachElevator;
     const char* _failureMsg;
 
-    // hard limits for cg from gear positions
+    float _cgMax;          // hard limits for cg from gear position
-    float _cgMaxX;
+    float _cgMin;          // hard limits for cg from gear position
-    float _cgMinX;
+    float _cgDesiredMax;   // desired cg max in %MAC from config
+    float _cgDesiredMin;   // desired cg min in %MAC from config
+    float _cgDesiredFront; // calculated desired cg x max
+    float _cgDesiredAft;   // calculated desired cg x min 
+    bool _autoBallast = false;
 };
 
 }; // namespace yasim

src/FDM/YASim/FGFDM.cpp

@@ -133,11 +133,12 @@ void FGFDM::init()
 
     // write some compile time information to property tree
     _yasimN->getNode("config-version",true)->setIntValue(_airplane.getVersion());
-    _yasimN->getNode("model/cg-x-min",true)->setFloatValue(_airplane.getCGMinX());
+    _yasimN->getNode("model/cg-x-min",true)->setFloatValue(_airplane.getCGHardLimitXMin());
-    _yasimN->getNode("model/cg-x-max",true)->setFloatValue(_airplane.getCGMaxX());
+    _yasimN->getNode("model/cg-x-max",true)->setFloatValue(_airplane.getCGHardLimitXMax());
     
     // prepare nodes for write at runtime
     _cg_x = _yasimN->getNode("cg-x-m", true);
+    _cg_xmacN = _yasimN->getNode("cg-x-mac", true);
     _cg_y = _yasimN->getNode("cg-y-m", true);
     _cg_z = _yasimN->getNode("cg-z-m", true);
     _vxN = _yasimN->getNode("velocities/v-x", true);
@@ -258,6 +259,10 @@ void FGFDM::startElement(const char* name, const XMLAttributes &atts)
       if( !_airplane.isVersionOrNewer( Version::YASIM_VERSION_CURRENT ) ) {
         SG_LOG(SG_FLIGHT, SG_DEV_ALERT, "This aircraft does not use the latest yasim configuration version.");
       }
+      _airplane.setDesiredCGRangeInPercentOfMAC(attrf(a, "cg-min", 0.1f), attrf(a, "cg-max", 0.3f)); //FIXME find reasonable defaults
+      if (attrb(a, "auto-ballast")) {
+        _airplane.setAutoBallast(true);
+      }
     } else if(eq(name, "approach")) {
       float spd, alt = 0;
       if (a->hasAttribute("speed")) { spd = attrf(a, "speed") * KTS2MPS; }
@@ -691,6 +696,7 @@ void FGFDM::setOutputProperties(float dt)
     _cg_x->setFloatValue(cg[0]);
     _cg_y->setFloatValue(cg[1]);
     _cg_z->setFloatValue(cg[2]);
+    _cg_xmacN->setFloatValue(_airplane.getCGMAC());
 
     State* s = _airplane.getModel()->getState();
     float v[3], acc[3], rot[3], racc[3];

src/FDM/YASim/FGFDM.hpp

@@ -118,6 +118,7 @@ private:
     SGPropertyNode_ptr _arxN;
     SGPropertyNode_ptr _aryN;
     SGPropertyNode_ptr _arzN;
+    SGPropertyNode_ptr _cg_xmacN;
 };
 
 }; // namespace yasim

src/FDM/YASim/Wing.cpp

@@ -213,6 +213,13 @@ void Wing::compile()
     _wingspan = Math::abs(2*_tip[1]);
     _aspectRatio = _wingspan / _meanChord;
 
+    _netSpan = Math::abs(2*(_tip[1]-_base[1]));
+    // http://www.nasascale.org/p2/wp-content/uploads/mac-calculator.htm
+    const float commonFactor = _chord*(0.5+_taper)/(3*_chord*(1+_taper));
+    _mac = _chord-(2*_chord*(1-_taper)*commonFactor);
+    _macRootDistance = _netSpan*commonFactor;
+    _macX = _base[0]-Math::tan(_sweep) * _macRootDistance + _mac/2;
+    
     // The wing's Y axis will be the "left" vector.  The Z axis will
     // be perpendicular to this and the local (!) X axis, because we
     // want motion along the local X axis to be zero AoA (i.e. in the

src/FDM/YASim/Wing.hpp

@@ -17,28 +17,28 @@ public:
 
     // Do we mirror ourselves about the XZ plane?
     void setMirror(bool mirror) { _mirror = mirror; }
-    bool isMirrored() { return _mirror; };
+    const bool isMirrored() { return _mirror; };
     
     // Wing geometry in local coordinates:
     
     // base point of wing
     void setBase(const float* base) { Math::set3(base, _base); }
-    void getBase(float* base) { Math::set3(_base, base); };
+    const void getBase(float* base) { Math::set3(_base, base); };
     // dist. ALONG wing (not span!)     
     void setLength(float length) { _length = length; }
-    float getLength() { return _length; };
+    const float getLength() { return _length; };
     // at base, measured along X axis
     void setChord(float chord) { _chord = chord; }
-    float getChord() { return _chord; };
+    const float getChord() { return _chord; };
     // fraction of chord at wing tip, 0..1
     void setTaper(float taper) { _taper = taper; }
-    float getTaper() { return _taper; };
+    const float getTaper() { return _taper; };
     // radians
     void setSweep(float sweep) { _sweep = sweep; }
-    float getSweep() { return _sweep; };
+    const float getSweep() { return _sweep; };
     // radians, positive is "up"
     void setDihedral(float dihedral) { _dihedral = dihedral; }
-    float getDihedral() { return _dihedral; };
+    const float getDihedral() { return _dihedral; };
     
     void setIncidence(float incidence);
     void setTwist(float angle) { _twist = angle; }
@@ -67,27 +67,31 @@ public:
 
     // Compile the thing into a bunch of Surface objects
     void compile();
-    void getTip(float* tip) { Math::set3(_tip, tip);};
+    const void getTip(float* tip) { Math::set3(_tip, tip);};
     
     // valid only after Wing::compile() was called
-    float getSpan() { return _wingspan; };
+    const float getSpan() { return _wingspan; };
-    float getArea() { return _wingspan*_meanChord; };
+    const float getArea() { return _wingspan*_meanChord; };
-    float getAspectRatio() { return _aspectRatio; };
+    const float getAspectRatio() { return _aspectRatio; };
-    float getSMC() { return _meanChord; };
+    const float getSMC() { return _meanChord; };
+    const float getMAC() { return _mac; }; // get length of MAC
+    const float getMACx() { return _macX; }; // get x-coord of MAC leading edge 
+    const float getMACy() { return _base[1]+_macRootDistance; }; // get y-coord of MAC leading edge 
     
-    int numSurfaces() { return _surfs.size(); }
+    
+    const int numSurfaces() { return _surfs.size(); }
     Surface* getSurface(int n) { return ((SurfRec*)_surfs.get(n))->surface; }
-    float getSurfaceWeight(int n) { return ((SurfRec*)_surfs.get(n))->weight; }
+    const float getSurfaceWeight(int n) { return ((SurfRec*)_surfs.get(n))->weight; }
 
     // The overall drag coefficient for the wing as a whole.  Units are
     // arbitrary.
     void setDragScale(float scale);
-    float getDragScale() { return _dragScale; }
+    const float getDragScale() { return _dragScale; }
 
     // The ratio of force along the Z (lift) direction of each wing
     // segment to that along the X (drag) direction.
     void setLiftRatio(float ratio);
-    float getLiftRatio() { return _liftRatio; }
+    const float getLiftRatio() { return _liftRatio; }
 
 private:
     void interp(const float* v1, const float* v2, const float frac, float* out);
@@ -114,6 +118,10 @@ private:
     // 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;
 

src/FDM/YASim/yasim-test.cpp

@@ -237,7 +237,9 @@ int main(int argc, char** argv)
     }
   }
   else {
-    printf("Solution results:");
+    printf("==========================\n");
+    printf("= YASim solution results =\n");
+    printf("==========================\n");
     float aoa = a->getCruiseAoA() * RAD2DEG;
     float tail = -1 * a->getTailIncidence() * RAD2DEG;
     float drag = 1000 * a->getDragCoefficient();
@@ -247,18 +249,26 @@ int main(int argc, char** argv)
 
     float SI_inertia[9];
     a->getModel()->getBody()->getInertiaMatrix(SI_inertia);
+    float MAC = a->getWing()->getMAC();
+    float MACx = a->getWing()->getMACx();
+    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\n", aoa);
+    printf("       Cruise AoA: %f deg\n", aoa);
-    printf("   Tail Incidence: %f\n", tail);
+    printf("   Tail Incidence: %f deg\n", tail);
-    printf("Approach Elevator: %f\n", a->getApproachElevator());
+    printf("Approach Elevator: %f\n\n", a->getApproachElevator());
-    printf("               CG: x:%.3f, y:%.3f, z:%.3f\n\n", cg[0], cg[1], cg[2]);
+    printf("               CG: x:%.3f, y:%.3f, z:%.3f\n", cg[0], cg[1], cg[2]);
-    printf("Inertia tensor [kg*m^2], origo at CG:\n");
+    printf("    Wing MAC (*1): x:%.2f, y:%.2f, length:%.1f \n", MACx, MACy, MAC);
+    printf("    CG-x rel. MAC: %.3f\n", a->getCGMAC());
+    printf("    CG-x  desired: %.3f < %.3f < %.3f \n", a->getCGSoftLimitXMin(), cg[0], a->getCGSoftLimitXMax());
+    
+    printf("\nInertia tensor [kg*m^2], origo at CG:\n\n");
     printf("  %7.3f, %7.3f, %7.3f\n", SI_inertia[0], SI_inertia[1], SI_inertia[2]);
     printf("  %7.3f, %7.3f, %7.3f\n", SI_inertia[3], SI_inertia[4], SI_inertia[5]);
     printf("  %7.3f, %7.3f, %7.3f\n", SI_inertia[6], SI_inertia[7], SI_inertia[8]);
+    printf("\n(*1) MAC calculation works on <wing> only! Numbers will be wrong for segmented wings, e.g. <wing>+<mstab>.\n");
   }
     delete fdm;
     return 0;