diff --git a/src/FDM/YASim/Airplane.cpp b/src/FDM/YASim/Airplane.cpp
index 2b32968f0..928d13448 100644
--- a/src/FDM/YASim/Airplane.cpp
+++ b/src/FDM/YASim/Airplane.cpp
@@ -16,6 +16,9 @@
 
 namespace yasim {
 
+//default prop names
+static const char* DEF_PROP_ELEVATOR_TRIM = "/controls/flight/elevator-trim";
+
 // gadgets
 inline float abs(float f) { return f<0 ? -f : f; }
 
@@ -59,12 +62,13 @@ Airplane::~Airplane()
     }
     for(i=0; i<_solveWeights.size(); i++)
         delete (SolveWeight*)_solveWeights.get(i);
-    for(i=0; i<_cruiseConfig.controls.size(); i++)
-        delete (ControlSetting*)_cruiseConfig.controls.get(i);
+    for(i=0; i<_cruiseConfig.controls.size(); i++) {
+        ControlSetting* c = (ControlSetting*)_cruiseConfig.controls.get(i);
+        delete c;        
+    }
     for(i=0; i<_approachConfig.controls.size(); i++) {
         ControlSetting* c = (ControlSetting*)_approachConfig.controls.get(i);
-        if(c != &_approachElevator)
-            delete c;
+        delete c;
     }
     delete _wing;
     delete _tail;
@@ -78,7 +82,6 @@ void Airplane::iterate(float dt)
 {
     // The gear might have moved.  Change their aerodynamics.
     updateGearState();
-
     _model.iterate();
 }
 
@@ -143,7 +146,8 @@ void Airplane::setApproach(float speed, float altitude, float aoa, float fuel, f
 {
     _approachConfig.speed = speed;
     _approachConfig.altitude = altitude;
-    _approachConfig.state.setupOrientationFromAoa(aoa); // see runConfig()
+    // solver assumes fixed (given) AoA for approach, so setup once
+    _approachConfig.state.setupOrientationFromAoa(aoa); 
     _approachConfig.aoa = aoa; // not strictly needed, see runConfig()
     _approachConfig.fuel = fuel;
     _approachConfig.glideAngle = gla;
@@ -153,23 +157,32 @@ void Airplane::setCruise(float speed, float altitude, float fuel, float gla)
 {
     _cruiseConfig.speed = speed;
     _cruiseConfig.altitude = altitude;
-    _cruiseConfig.aoa = 0;
-    _tailIncidence.val = 0;
     _cruiseConfig.fuel = fuel;
     _cruiseConfig.glideAngle = gla;
 }
 
-void Airplane::setElevatorControl(const char* prop)
+/// set property name for elevator  
+void Airplane::setElevatorControl(const char* propName)
 {
-    _approachElevator.propHandle = getControlMap()->getInputPropertyHandle(prop);
-    _approachElevator.val = 0;
-    _approachConfig.controls.add(&_approachElevator);
+    _approachElevator = _addControlSetting(APPROACH, propName, 0);
+}
+
+/// set property name for hstab trim
+void Airplane::setHstabTrimControl(const char* propName)
+{
+    _tailIncidence = _addControlSetting(APPROACH, propName, 0);
+    _ti2 = _addControlSetting(CRUISE, propName, 0);
 }
 
 void Airplane::addControlSetting(Configuration cfg, const char* prop, float val)
+{
+    _addControlSetting(cfg, prop,val);
+}
+
+Airplane::ControlSetting* Airplane::_addControlSetting(Configuration cfg, const char* prop, float val)
 {
     ControlSetting* c = new ControlSetting();
-    c->propHandle = getControlMap()->getInputPropertyHandle(prop);
+    c->propHandle = _controlMap.getInputPropertyHandle(prop);
     c->val = val;
     switch (cfg) {
         case APPROACH:
@@ -179,6 +192,21 @@ void Airplane::addControlSetting(Configuration cfg, const char* prop, float val)
             _cruiseConfig.controls.add(c);
             break;
     }
+    return c;
+}
+
+
+/**
+ * used by the XML parser in FGFDM and solveAirplane
+ */
+void Airplane::addControlInput(const char* propName, ControlMap::ControlType type, void* obj, int subobj, int opt, float src0, float src1, float dst0, float dst1)
+{
+    ControlMap::ObjectID oid = ControlMap::getObjectID(obj, subobj);    
+    _controlMap.addMapping(propName, type, oid, opt, src0, src1, dst0, dst1);
+    // tail incidence is needed by solver so capture the prop name if used in XML
+    if (type == ControlMap::INCIDENCE && obj == _tail) {
+        setHstabTrimControl(propName);
+    }
 }
 
 void Airplane::addSolutionWeight(Configuration cfg, int idx, float wgt)
@@ -727,12 +755,14 @@ void Airplane::setupWeights(bool isApproach)
     }
 }
 
+/// used by solver to simulate property input
 void Airplane::setControlValues(const Vector& controls)
 {
     _controlMap.reset();
     for(int i=0; i < controls.size(); i++) {
         ControlSetting* c = (ControlSetting*)controls.get(i);
-        _controlMap.setInput(c->propHandle, c->val);
+        if (c->propHandle >= 0)
+            _controlMap.setInput(c->propHandle, c->val);
     }
     _controlMap.applyControls(); 
 }
@@ -740,7 +770,7 @@ void Airplane::setControlValues(const Vector& controls)
 void Airplane::runConfig(Config &cfg)
 {
     // aoa is consider to be given for approach so we calculate orientation 
-    // only once in setApproach()
+    // for approach only once in setApproach() but everytime for cruise here.
     if (!cfg.isApproach) {
         cfg.state.setupOrientationFromAoa(cfg.aoa);
     }
@@ -774,6 +804,7 @@ void Airplane::runConfig(Config &cfg)
     _model.initIteration();
     _model.calcForces(&cfg.state);
 }
+
 /// Used only in solveAirplane() and solveHelicopter(), not at runtime
 void Airplane::applyDragFactor(float factor)
 {
@@ -790,23 +821,22 @@ void Airplane::applyDragFactor(float factor)
     }
     for(i=0; i<_fuselages.size(); i++) {
         Fuselage* f = (Fuselage*)_fuselages.get(i);
-        int j;
-        for(j=0; j<f->surfs.size(); j++) {
+        for(int j=0; j<f->surfs.size(); j++) {
             Surface* s = (Surface*)f->surfs.get(j);
             if( isVersionOrNewer( YASIM_VERSION_32 ) ) {
-            // For new YASim, the solver drag factor is only applied to
-            // the X axis for Fuselage Surfaces.
-            // The solver is tuning the coefficient for longitudinal drag,
-            // along the direction of flight. A fuselage's lateral drag
-            // is completely independent and is normally much higher;
-            // it won't be affected by the streamlining done to reduce
-            // longitudinal drag. So the solver should only adjust the
-            // fuselage's longitudinal (X axis) drag coefficient.
-            s->setDragCoefficient(s->getDragCoefficient() * applied);
+                // For new YASim, the solver drag factor is only applied to
+                // the X axis for Fuselage Surfaces.
+                // The solver is tuning the coefficient for longitudinal drag,
+                // along the direction of flight. A fuselage's lateral drag
+                // is completely independent and is normally much higher;
+                // it won't be affected by the streamlining done to reduce
+                // longitudinal drag. So the solver should only adjust the
+                // fuselage's longitudinal (X axis) drag coefficient.
+                s->mulDragCoefficient(applied);
             } else {
-            // Originally YASim applied the drag factor to all axes
-            // for Fuselage Surfaces.
-            s->mulTotalForceCoefficient(applied);
+                // Originally YASim applied the drag factor to all axes
+                // for Fuselage Surfaces.
+                s->mulTotalForceCoefficient(applied);
             }
         }
     }
@@ -820,7 +850,8 @@ void Airplane::applyDragFactor(float factor)
     }
 }
 
-/// Used only in Airplane::solve() and solveHelicopter(), not at runtime
+/// Used only in solveAirplane() and solveHelicopter(), not at runtime
+/// change lift coefficient cz in surfaces
 void Airplane::applyLiftRatio(float factor)
 {
     float applied = Math::pow(factor, SOLVE_TWEAK);
@@ -870,9 +901,20 @@ void Airplane::solveAirplane()
     _solutionIterations = 0;
     _failureMsg = 0;
 
+    if (_approachElevator == nullptr) {        
+        setElevatorControl(DEF_PROP_ELEVATOR_TRIM);
+    }
+    
+    if (_tailIncidence == nullptr) {
+        // no control mapping from XML parser, so we just create "local" 
+        // variables for solver instead of full mapping / property
+        _tailIncidence = new ControlSetting;
+        _ti2 = new ControlSetting;
+    }
+
     while(1) {
-        if(_solutionIterations++ > 10000) { 
-            _failureMsg = "Solution failed to converge after 10000 iterations";
+        if(_solutionIterations++ > SOLVER_MAX_ITERATIONS) { 
+            _failureMsg = "Solution failed to converge!";
             return;
         }
         // Run an iteration at cruise, and extract the needed numbers:
@@ -884,7 +926,7 @@ void Airplane::solveAirplane()
         _cruiseConfig.state.localToGlobal(tmp, tmp);
         float xforce = _cruiseConfig.weight * tmp[0];
         float clift0 = _getLift(_cruiseConfig);
-        float pitch0 = _getPitch(_cruiseConfig);
+        float cpitch0 = _getPitch(_cruiseConfig);
 
         // Run an approach iteration, and do likewise
         runConfig(_approachConfig);
@@ -899,19 +941,23 @@ void Airplane::solveAirplane()
         float clift1 = _getLift(_cruiseConfig);
 
         // Do the same with the tail incidence
-        _tail->setIncidence(_tailIncidence.val + ARCMIN);
+        float savedIncidence = _tailIncidence->val;
+        _ti2->val = _tailIncidence->val += ARCMIN;
+        if (!_tail->setIncidence(_tailIncidence->val)) {
+            _failureMsg = "Tail incidence out of bounds.";
+            return;
+        };
         runConfig(_cruiseConfig);
-        _tail->setIncidence(_tailIncidence.val);
+        _ti2->val = _tailIncidence->val = savedIncidence;
+        _tail->setIncidence(_tailIncidence->val);
 
-        float pitch1 = _getPitch(_cruiseConfig);
+        float cpitch1 = _getPitch(_cruiseConfig);
 
         // Now calculate:
         float awgt = 9.8f * _approachConfig.weight;
 
         float dragFactor = thrust / (thrust-xforce);
         float liftFactor = awgt / (awgt+alift);
-        float aoaDelta = -clift0 * (ARCMIN/(clift1-clift0));
-        float tailDelta = -pitch0 * (ARCMIN/(pitch1-pitch0));
 
         // Sanity:
         if(dragFactor <= 0 || liftFactor <= 0)
@@ -922,12 +968,11 @@ void Airplane::solveAirplane()
         // same thing -- pitching moment -- by diddling a different
         // variable).
         const float ELEVDIDDLE = 0.001f;
-        _approachElevator.val += ELEVDIDDLE;
+        _approachElevator->val += ELEVDIDDLE;
         runConfig(_approachConfig);
-        _approachElevator.val -= ELEVDIDDLE;
+        _approachElevator->val -= ELEVDIDDLE;
 
         double apitch1 = _getPitch(_approachConfig);
-        float elevDelta = -apitch0 * (ELEVDIDDLE/(apitch1-apitch0));
 
         // Now apply the values we just computed.  Note that the
         // "minor" variables are deferred until we get the lift/drag
@@ -944,25 +989,28 @@ void Airplane::solveAirplane()
         }
 
         // OK, now we can adjust the minor variables:
+        float aoaDelta = -clift0 * (ARCMIN/(clift1-clift0));
+        float tailDelta = -cpitch0 * (ARCMIN/(cpitch1-cpitch0));
         _cruiseConfig.aoa += SOLVE_TWEAK*aoaDelta;
-        _tailIncidence.val += SOLVE_TWEAK*tailDelta;
+        _tailIncidence->val += SOLVE_TWEAK*tailDelta;
         
         _cruiseConfig.aoa = Math::clamp(_cruiseConfig.aoa, -0.175f, 0.175f);
-        _tailIncidence.val = Math::clamp(_tailIncidence.val, -0.175f, 0.175f);
+        _tailIncidence->val = Math::clamp(_tailIncidence->val, -0.175f, 0.175f);
 
         if(abs(xforce/_cruiseConfig.weight) < STHRESH*0.0001 &&
-        abs(alift/_approachConfig.weight) < STHRESH*0.0001 &&
-        abs(aoaDelta) < STHRESH*.000017 &&
-        abs(tailDelta) < STHRESH*.000017)
+          abs(alift/_approachConfig.weight) < STHRESH*0.0001 &&
+          abs(aoaDelta) < STHRESH*.000017 &&
+          abs(tailDelta) < STHRESH*.000017)
         {
+            float elevDelta = -apitch0 * (ELEVDIDDLE/(apitch1-apitch0));
             // If this finaly value is OK, then we're all done
             if(abs(elevDelta) < STHRESH*0.0001)
                 break;
 
             // Otherwise, adjust and do the next iteration
-            _approachElevator.val += SOLVE_TWEAK * elevDelta;
-            if(abs(_approachElevator.val) > 1) {
-                _failureMsg = "Insufficient elevator to trim for approach";
+            _approachElevator->val += SOLVE_TWEAK * elevDelta;
+            if(abs(_approachElevator->val) > 1) {
+                _failureMsg = "Insufficient elevator to trim for approach.";
                 break;
             }
         }
@@ -977,10 +1025,16 @@ void Airplane::solveAirplane()
     } else if(Math::abs(_cruiseConfig.aoa) >= .17453293) {
         _failureMsg = "Cruise AoA > 10 degrees";
         return;
-    } else if(Math::abs(_tailIncidence.val) >= .17453293) {
+    } else if(Math::abs(_tailIncidence->val) >= .17453293) {
         _failureMsg = "Tail incidence > 10 degrees";
         return;
     }
+    // if we have a property tree, export result from solver
+    if (_wingsN != nullptr) {
+        if (_tailIncidence->propHandle >= 0) {
+            fgSetFloat(_controlMap.getProperty(_tailIncidence->propHandle)->name, _tailIncidence->val);
+        }
+    }
 }
 
 void Airplane::solveHelicopter()
diff --git a/src/FDM/YASim/Airplane.hpp b/src/FDM/YASim/Airplane.hpp
index 3fa34d471..5882a4a15 100644
--- a/src/FDM/YASim/Airplane.hpp
+++ b/src/FDM/YASim/Airplane.hpp
@@ -63,11 +63,11 @@ public:
 
     void setApproach(float speed, float altitude, float aoa, float fuel, float gla);
     void setCruise(float speed, float altitude, float fuel, float gla);
-
-     ///set name of property controlling the elevator
-    void setElevatorControl(const char* prop);
+    
+    /// add (fixed) control setting to approach/cruise config (for solver)
     void addControlSetting(Configuration cfg, const char* prop, float val);
-
+    /// add a control input mapping for runtime
+    void addControlInput(const char* propName, ControlMap::ControlType type, void* obj, int subobj, int opt, float src0, float src1, float dst0, float dst1);
     void addSolutionWeight(Configuration cfg, int idx, float wgt);
 
     int numGear() const { return _gears.size(); }
@@ -101,10 +101,11 @@ public:
     float getDragCoefficient() const { return _dragFactor; }
     float getLiftRatio() const { return _liftRatio; }
     float getCruiseAoA() const { return _cruiseConfig.aoa; }
-    float getTailIncidence() const { return _tailIncidence.val; }
-    float getApproachElevator() const { return _approachElevator.val; }
+    float getTailIncidence() const { return _tailIncidence->val; }
+    float getApproachElevator() const { return _approachElevator->val; }
     const char* getFailureMsg() const { return _failureMsg; }
 
+    // next two are used only in yasim CLI tool
     void setApproachControls() { setControlValues(_approachConfig.controls); }
     void setCruiseControls() { setControlValues(_cruiseConfig.controls); }
     
@@ -208,6 +209,11 @@ private:
     float _getWingLoad(float mass) const;
     ///calculate distance between CGx and AC of wing w
     float _getWingLever(const Wing* w) const;
+    ControlSetting* _addControlSetting(Configuration cfg, const char* prop, float val);
+     ///set name of property controlling the elevator
+    void setElevatorControl(const char* propName);
+    /// set property name controling tail trim (incidence)
+    void setHstabTrimControl(const char* propName);
     
     Model _model;
     ControlMap _controlMap;
@@ -237,8 +243,9 @@ private:
     int _solutionIterations {0};
     float _dragFactor {1};
     float _liftRatio {1};
-    ControlSetting _tailIncidence;
-    ControlSetting _approachElevator;
+    ControlSetting* _tailIncidence {nullptr}; // added to approach config so solver can change it
+    ControlSetting* _ti2 {nullptr}; // copy of _tailIncidence added to cruise config
+    ControlSetting* _approachElevator {nullptr};
     const char* _failureMsg {0};
     /// hard limits for cg from gear position
     float _cgMax {-1e6};         
diff --git a/src/FDM/YASim/ControlMap.cpp b/src/FDM/YASim/ControlMap.cpp
index f1aac7ef8..15d3ce9c2 100644
--- a/src/FDM/YASim/ControlMap.cpp
+++ b/src/FDM/YASim/ControlMap.cpp
@@ -88,10 +88,11 @@ ControlMap::~ControlMap()
 }
 
 /**
-prop: name of input property
+inputProp: name of input property
 control: identifier (see enum OutputType)
-object: object to which this input belongs to
+id: object to which this input belongs to
 options: bits OPT_INVERT, OPT_SPLIT, OPT_SQUARE
+src,dst: input will be clamped to src range and mapped to dst range 
 */
 void ControlMap::addMapping(const char* inputProp, ControlType control, ObjectID id, int options, float src0, float src1, float dst0, float dst1)
 {
@@ -150,9 +151,9 @@ void ControlMap::reset()
     }
 }
 
-void ControlMap::setInput(int input, float val)
+void ControlMap::setInput(int propHandle, float val)
 {
-    Vector* maps = (Vector*)_inputs.get(input);
+    Vector* maps = (Vector*)_inputs.get(propHandle);
     for(int i = 0; i < maps->size(); i++) {
         MapRec* m = (MapRec*)maps->get(i);
         float val2 = val;
@@ -186,28 +187,22 @@ float ControlMap::getOutputR(int handle)
 
 void ControlMap::applyControls(float dt)
 {
-    int outrec;
-    for(outrec=0; outrec<_outputs.size(); outrec++) 
+    for(int outrec=0; outrec<_outputs.size(); outrec++) 
     {
         OutRec* o = (OutRec*)_outputs.get(outrec);
     
         // Generate a summed value.  Note the check for "split"
         // control axes like ailerons.
         float lval = 0, rval = 0;
-        int i;
-        for(i=0; i<o->maps.size(); i++) {
+        for(int i = 0; i < o->maps.size(); i++) {
             MapRec* m = (MapRec*)o->maps.get(i);
             float val = m->val;
 
-            if(m->opt & OPT_SQUARE)
-            val = val * Math::abs(val);
-            if(m->opt & OPT_INVERT)
-            val = -val;
+            if(m->opt & OPT_SQUARE) { val = val * Math::abs(val); }
+            if(m->opt & OPT_INVERT) { val = -val; }
             lval += val;
-            if(m->opt & OPT_SPLIT)
-            rval -= val;
-            else
-            rval += val;
+            if(m->opt & OPT_SPLIT) { rval -= val; }
+            else { rval += val; }
         }
 
         // If there is a finite transition time, clamp the values to
@@ -218,12 +213,15 @@ void ControlMap::applyControls(float dt)
             float adl = Math::abs(dl);
             float adr = Math::abs(dr);
         
-            float max = (dt/o->transitionTime) * (rangeMax(o->control) - rangeMin(o->control));
-            if(adl > max) dl = dl*max/adl;
-            if(adr > max) dr = dr*max/adr;
-
-            lval = o->oldValueLeft + dl;
-            rval = o->oldValueRight + dr;
+            float maxDelta = (dt/o->transitionTime) * (rangeMax(o->control) - rangeMin(o->control));
+            if(adl > maxDelta) {
+                dl = dl*maxDelta/adl;
+                lval = o->oldValueLeft + dl;
+            }
+            if(adr > maxDelta) {
+                dr = dr*maxDelta/adr;
+                rval = o->oldValueRight + dr;
+            }
         }
 
         o->oldValueLeft = lval;
@@ -357,6 +355,7 @@ void ControlMap::applyControls(float dt)
             case PROP:
                 break;
             case INCIDENCE:
+                ((Wing*)obj)->setIncidence(lval);
                 break;
         }
     }
@@ -366,6 +365,7 @@ float ControlMap::rangeMin(ControlType control)
 {
     // The minimum of the range for each type of control
     switch(control) {
+        case INCIDENCE: return INCIDENCE_MIN;
         case FLAP0:    return -1;  // [-1:1]
         case FLAP1:    return -1;
         case STEER:    return -1;
@@ -384,6 +384,7 @@ float ControlMap::rangeMax(ControlType control)
 {
     // The maximum of the range for each type of control
     switch(control) {
+        case INCIDENCE: return INCIDENCE_MAX;
         case FLAP0:    return 1; // [-1:1]
         case FLAP1:    return 1;
         case STEER:    return 1;
@@ -397,13 +398,14 @@ float ControlMap::rangeMax(ControlType control)
 /// register property name, return ID (int)
 int ControlMap::getInputPropertyHandle(const char* name)
 {
+    // search for existing
     for(int i=0; i < _properties.size(); i++) {
         PropHandle* p = (PropHandle*)_properties.get(i);
         if(!strcmp(p->name, name))
             return p->handle;
     }
 
-    // create new
+    // else create new
     PropHandle* p = new PropHandle();
     p->name = strdup(name);
     
@@ -411,6 +413,7 @@ int ControlMap::getInputPropertyHandle(const char* name)
 
     Vector* v = new Vector();
     p->handle = _inputs.add(v);
+
     _properties.add(p);
     return p->handle;
 }
@@ -446,4 +449,10 @@ ControlMap::ObjectID ControlMap::getObjectID(void* object, int subObj)
     return o;
 }
 
+// used at runtime in FGFDM::getExternalInput
+ControlMap::PropHandle* ControlMap::getProperty(const int i) {
+    assert((i >= 0) && (i < _properties.size()));
+    return ((PropHandle*)_properties.get(i));
+}
+
 } // namespace yasim
diff --git a/src/FDM/YASim/ControlMap.hpp b/src/FDM/YASim/ControlMap.hpp
index 9ebecd22e..68916a1f3 100644
--- a/src/FDM/YASim/ControlMap.hpp
+++ b/src/FDM/YASim/ControlMap.hpp
@@ -72,6 +72,7 @@ public:
         char* name {nullptr};
         int handle {0};
     };
+    // to identify controls per wing section we need wing object + section id 
     struct ObjectID {
         void* object {nullptr};
         int subObj {0};
@@ -125,7 +126,7 @@ public:
     // register property name, return handle
     int getInputPropertyHandle(const char* name);
     int numProperties() { return _properties.size(); }
-    PropHandle* getProperty(const int i) { return ((PropHandle*)_properties.get(i)); }
+    PropHandle* getProperty(const int i);
 
 private:
     //output data for a control of an object
@@ -154,10 +155,9 @@ private:
 
     // An unordered list of output settings.
     Vector _outputs;
-    // control properties
-    Vector _properties;
+    
+    Vector _properties; // list of PropHandle*
 
-    void* addMapping(const char* prop, ControlType control, ObjectID id, int options = 0);
     OutRec* getOutRec(ObjectID id, ControlType control);
 };
 
diff --git a/src/FDM/YASim/FGFDM.cpp b/src/FDM/YASim/FGFDM.cpp
index 043cdec06..a5ab68121 100644
--- a/src/FDM/YASim/FGFDM.cpp
+++ b/src/FDM/YASim/FGFDM.cpp
@@ -37,11 +37,6 @@ namespace yasim {
 
 FGFDM::FGFDM()
 {
-    // Map /controls/flight/elevator to the approach elevator control.  This
-    // should probably be settable, but there are very few aircraft
-    // who trim their approaches using things other than elevator.
-    _airplane.setElevatorControl("/controls/flight/elevator-trim");
-
     // FIXME: read seed from somewhere?
     int seed = 0;
     _turb = new Turbulence(10, seed);
@@ -62,8 +57,8 @@ FGFDM::~FGFDM()
         delete wr;
     }
 
-    for(int i=0; i<_controlProps.size(); i++)
-        delete (PropOut*)_controlProps.get(i);
+    for(int i=0; i<_controlOutputs.size(); i++)
+        delete (ControlOutput*)_controlOutputs.get(i);
 
     delete _turb;
 }
@@ -419,14 +414,15 @@ void FGFDM::setOutputProperties(float dt)
     _arzN->setFloatValue(racc[2]);
 
     ControlMap* cm = _airplane.getControlMap();
-    for(int i=0; i<_controlProps.size(); i++) {
-        PropOut* p = (PropOut*)_controlProps.get(i);
+    for(int i=0; i<_controlOutputs.size(); i++) {
+        ControlOutput* p = (ControlOutput*)_controlOutputs.get(i);
         float val = (p->left
                      ? cm->getOutput(p->handle)
                      : cm->getOutputR(p->handle));
         float rmin = cm->rangeMin(p->control);
         float rmax = cm->rangeMax(p->control);
         float frac = (val - rmin) / (rmax - rmin);
+        // clamp output
         val = frac*(p->max - p->min) + p->min;
         p->prop->setFloatValue(val);
     }
@@ -531,11 +527,10 @@ void FGFDM::parseWing(const XMLAttributes* a, const char* type, Airplane* airpla
     float chord {0};
     float length = attrf(a, "length");
 
-    // These come in with positive indicating positive AoA, but the
-    // internals expect a rotation about the left-pointing Y axis, so
-    // invert the sign.
+    // positive incidence/twist means positive AoA (leading edge up). 
+    // Due to the coordinate system used in class Surface the sign will be inverted (only) there.
     float incidence {0};
-    float twist = attrf(a, "twist", 0) * DEG2RAD * -1;
+    float twist = attrf(a, "twist", 0) * DEG2RAD;
     
     // if this element is declared as section of a wing, skip attributes 
     // that are ignored in class Wing anyway because they are calculated
@@ -543,7 +538,15 @@ void FGFDM::parseWing(const XMLAttributes* a, const char* type, Airplane* airpla
     if (!isSection) {
         attrf_xyz(a, base);
         chord = attrf(a, "chord");
-        incidence = attrf(a, "incidence", 0) * DEG2RAD * -1;
+        incidence = attrf(a, "incidence", 0) * DEG2RAD;
+    }
+    else {
+        if (a->hasAttribute("x") || a->hasAttribute("y") || a->hasAttribute("z") ||
+            a->hasAttribute("chord") || a->hasAttribute("incidence") 
+        ) {
+            SG_LOG(SG_FLIGHT, SG_WARN, "YASim warning: redundant attribute in wing definition \n" 
+            "when using <wing append=\"1\" ...> x, y, z, chord and incidence will be ignored. ");
+        }
     } 
 
     // optional attributes (with defaults)
@@ -569,8 +572,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);
+        if (a->hasAttribute("incidence-max-deg")) w->setIncidenceMax(attrf(a, "incidence-max-deg") * DEG2RAD);
+        if (a->hasAttribute("incidence-min-deg")) w->setIncidenceMin(attrf(a, "incidence-min-deg") * DEG2RAD);
     } else {
         w = new Wing(airplane, mirror);
     }
@@ -1155,48 +1158,47 @@ void FGFDM::parseControlSetting(const XMLAttributes* a)
 void FGFDM::parseControlIn(const XMLAttributes* a)
 {
     // map input property to a YASim control 
-    ControlMap* cm = _airplane.getControlMap();
-    ControlMap::ControlType control = cm->parseControl(a->getValue("control"));
-    ControlMap::ObjectID oid = cm->getObjectID(_currObj, _wingSection);
+    ControlMap::ControlType control = ControlMap::parseControl(a->getValue("control"));
     int opt = 0;
     opt |= a->hasAttribute("split") ? ControlMap::OPT_SPLIT : 0;
     opt |= a->hasAttribute("invert") ? ControlMap::OPT_INVERT : 0;
     opt |= a->hasAttribute("square") ? ControlMap::OPT_SQUARE : 0;
     float src0, src1, dst0, dst1;
-    src0 = dst0 = cm->rangeMin(control);
-    src1 = dst1 = cm->rangeMax(control);
+    src0 = dst0 = ControlMap::rangeMin(control);
+    src1 = dst1 = ControlMap::rangeMax(control);
     if(a->hasAttribute("src0")) {
         src0 = attrf(a, "src0");
         src1 = attrf(a, "src1");
         dst0 = attrf(a, "dst0");
         dst1 = attrf(a, "dst1");
     }
-    cm->addMapping(a->getValue("axis"), control, oid, opt, src0, src1, dst0, dst1);
+    _airplane.addControlInput(a->getValue("axis"), control, _currObj, _wingSection, opt, src0, src1, dst0, dst1);
 }
 
 void FGFDM::parseControlOut(const XMLAttributes* a)
 {
     // A property output for a control on the current object
     ControlMap* cm = _airplane.getControlMap();
-    ControlMap::ControlType control = cm->parseControl(a->getValue("control"));
-    ControlMap::ObjectID oid = cm->getObjectID(_currObj, _wingSection);
+    ControlMap::ControlType control = ControlMap::parseControl(a->getValue("control"));
+    ControlMap::ObjectID oid = ControlMap::getObjectID(_currObj, _wingSection);
 
-    PropOut* p = new PropOut();
+    ControlOutput* p = new ControlOutput();
     p->prop = fgGetNode(a->getValue("prop"), true);
     p->handle = cm->getOutputHandle(oid, control);
     p->control = control;
     p->left = !(a->hasAttribute("side") &&
                     !strcmp("right", a->getValue("side")));
+    // for output clamping
     p->min = attrf(a, "min", cm->rangeMin(control));
     p->max = attrf(a, "max", cm->rangeMax(control));
-    _controlProps.add(p);
+    _controlOutputs.add(p);
 }
 
 void FGFDM::parseControlSpeed(const XMLAttributes* a)
 {
     ControlMap* cm = _airplane.getControlMap();
-    ControlMap::ControlType control = cm->parseControl(a->getValue("control"));
-    ControlMap::ObjectID oid = cm->getObjectID(_currObj, _wingSection);
+    ControlMap::ControlType control = ControlMap::parseControl(a->getValue("control"));
+    ControlMap::ObjectID oid = ControlMap::getObjectID(_currObj, _wingSection);
     int handle = cm->getOutputHandle(oid, control);
     float time = attrf(a, "transition-time", 0);
     cm->setTransitionTime(handle, time);
diff --git a/src/FDM/YASim/FGFDM.hpp b/src/FDM/YASim/FGFDM.hpp
index 6983910bf..18dbac131 100644
--- a/src/FDM/YASim/FGFDM.hpp
+++ b/src/FDM/YASim/FGFDM.hpp
@@ -40,9 +40,9 @@ private:
         float size {0};
         int handle {0};
     };
-    struct PropOut { 
-        SGPropertyNode* prop {nullptr};
+    struct ControlOutput { 
         int handle {0};
+        SGPropertyNode* prop {nullptr};
         ControlMap::ControlType control;
         bool left {false};
         float min {0};
@@ -107,7 +107,7 @@ private:
     Vector _thrusters;
 
     // Output properties for the ControlMap
-    Vector _controlProps;
+    Vector _controlOutputs;
 
     // Radius of the vehicle, for intersection testing.
     float _vehicle_radius {0};
diff --git a/src/FDM/YASim/Surface.cpp b/src/FDM/YASim/Surface.cpp
index f413f13ac..dc21656de 100644
--- a/src/FDM/YASim/Surface.cpp
+++ b/src/FDM/YASim/Surface.cpp
@@ -30,9 +30,9 @@ Surface::Surface(Version* version, const float* pos, float c0 = 1 ) :
         _flapN = _surfN->getNode("flap-pos", true);
         _slatN = _surfN->getNode("slat-pos", true);
         _spoilerN = _surfN->getNode("spoiler-pos", true);
-        _incidenceN = _surfN->getNode("incidence", true);
+        _incidenceN = _surfN->getNode("incidence-deg", true);
         _incidenceN->setFloatValue(0);
-        _twistN = _surfN->getNode("twist", true);
+        _twistN = _surfN->getNode("twist-deg", true);
         _twistN->setFloatValue(0);
         _surfN->getNode("pos-x", true)->setFloatValue(pos[0]);
         _surfN->getNode("pos-y", true)->setFloatValue(pos[1]);
@@ -340,17 +340,17 @@ float Surface::controlDrag(float lift, float drag)
 
 
 void Surface::setIncidence(float angle) {
-    _incidence = angle;
+    _incidence = angle * -1;
     if (_surfN != 0) {
-        _incidenceN->setFloatValue(angle);
+        _incidenceN->setFloatValue(angle * RAD2DEG);
     }
 }
 
 
 void Surface::setTwist(float angle) {
-    _twist = angle;
+    _twist = angle * -1;
     if (_surfN != 0) {
-        _twistN->setFloatValue(angle);
+        _twistN->setFloatValue(angle * RAD2DEG);
     }
 }
 
diff --git a/src/FDM/YASim/Wing.cpp b/src/FDM/YASim/Wing.cpp
index b837b2f17..aec0e8d12 100644
--- a/src/FDM/YASim/Wing.cpp
+++ b/src/FDM/YASim/Wing.cpp
@@ -390,19 +390,23 @@ void Wing::multiplyDragCoefficient(float factor)
 }
 
 ///update incidence for wing (rotate wing while maintaining initial twist config)
-void Wing::setIncidence(float incidence)
+bool Wing::setIncidence(float incidence)
 {
     if (incidence < _incidenceMin || incidence > _incidenceMax)
     {
-        fprintf(stderr, "YASim: cannot set incidence, parameter out of range.");
-        return;
+        return false;
     }
+    _incidence = incidence;
     WingSection* ws;
     for (int section=0; section < _sections.size(); section++) 
     {        
         ws = (WingSection*)_sections.get(section);
         ws->setIncidence(incidence);
     }    
+    if (_wingN) {
+        _wingN->getNode("incidence-deg", true)->setFloatValue(_incidence * RAD2DEG);
+    }
+    return true;
 }
     
 void Wing::WingSection::setDragCoefficient(float scale)
@@ -542,7 +546,7 @@ void Wing::writeInfoToProptree()
         sectN->getNode("base-y", true)->setFloatValue(ws->_rootChord.y);
         sectN->getNode("base-z", true)->setFloatValue(ws->_rootChord.z);
         sectN->getNode("chord", true)->setFloatValue(ws->_rootChord.length);
-        sectN->getNode("incidence", true)->setFloatValue(ws->_sectionIncidence);
+        sectN->getNode("incidence-deg", true)->setFloatValue(ws->_sectionIncidence * RAD2DEG);
     }
     _wingN->getNode("weight", true)->setFloatValue(wgt);
     _wingN->getNode("drag", true)->setFloatValue(dragSum);
@@ -596,8 +600,8 @@ void Wing::printSectionInfo()
     printf("#wing sections: %d\n", _sections.size());
     for (int section=0; section < _sections.size(); section++) {
         ws = (WingSection*)_sections.get(section);
-        printf("Section %d base point (%.3f, %.3f, %.3f), chord %.3f\n", section,
-               ws->_rootChord.x, ws->_rootChord.y, ws->_rootChord.z, ws->_rootChord.length);
+        printf("Section %d base point (%.3f, %.3f, %.3f), chord %.3f, incidence at section root %.1fdeg\n", section,
+               ws->_rootChord.x, ws->_rootChord.y, ws->_rootChord.z, ws->_rootChord.length, ws->_sectionIncidence * RAD2DEG);
     }
 }
 
diff --git a/src/FDM/YASim/Wing.hpp b/src/FDM/YASim/Wing.hpp
index c755b3b55..8b7c6bf71 100644
--- a/src/FDM/YASim/Wing.hpp
+++ b/src/FDM/YASim/Wing.hpp
@@ -164,7 +164,7 @@ public:
     void multiplyLiftRatio(float factor);
     void multiplyDragCoefficient(float factor);
     // setIncidence used to rotate (trim) the hstab
-    void setIncidence(float incidence);
+    bool setIncidence(float incidence);
     // limits for setIncidence
     void setIncidenceMin(float min) { _incidenceMin = min; };
     void setIncidenceMax(float max) { _incidenceMax = max; };
diff --git a/src/FDM/YASim/yasim-common.hpp b/src/FDM/YASim/yasim-common.hpp
index c591d19ba..0d0027776 100644
--- a/src/FDM/YASim/yasim-common.hpp
+++ b/src/FDM/YASim/yasim-common.hpp
@@ -8,6 +8,7 @@
  common file for YASim wide constants and static helper functions
  */
 namespace yasim {
+    static const int SOLVER_MAX_ITERATIONS = 4000;
     static const float YASIM_PI = 3.14159265358979323846f;
     static const float PI2 = YASIM_PI*2;
     static const float RAD2DEG = 180/YASIM_PI;
diff --git a/src/FDM/YASim/yasim-test.cpp b/src/FDM/YASim/yasim-test.cpp
index 1c6dc8276..53c70808d 100644
--- a/src/FDM/YASim/yasim-test.cpp
+++ b/src/FDM/YASim/yasim-test.cpp
@@ -36,7 +36,7 @@ enum Config
 };
 
 // Generate a graph of lift, drag and L/D against AoA at the specified
-// speed and altitude.  The result is a space-separated file of
+// speed and altitude.  The result is a tab-separated file of
 // numbers: "aoa lift drag LD" (aoa in degrees, lift and drag in
 // G's).  You can use this in gnuplot like so (assuming the output is
 // in a file named "dat":
@@ -95,11 +95,11 @@ void yasim_graph(Airplane* a, const float alt, const float kts, int cfg = CONFIG
       ld_max= ld;
       ld_max_deg = deg;
     }    
-    printf("%d %.4g %.4g %.4g\n", deg, lift, drag, ld);
+    printf("%2d\t%.4f\t%.4f\t%.4f\n", deg, lift, drag, ld);
   }
-  printf("# cl_max %g at %d deg\n", cl_max, cl_max_deg);
-  printf("# cd_min %g at %d deg\n", cd_min, cd_min_deg);
-  printf("# ld_max %g at %d deg\n", ld_max, ld_max_deg);  
+  printf("# cl_max %.4f at %d deg\n", cl_max, cl_max_deg);
+  printf("# cd_min %.4f at %d deg\n", cd_min, cd_min_deg);
+  printf("# ld_max %.4f at %d deg\n", ld_max, ld_max_deg);  
 }
 
 void yasim_masses(Airplane* a)
@@ -164,80 +164,13 @@ void yasim_drag(Airplane* a, const float aoa, const float alt, int cfg = CONFIG_
   printf("# cd_min %g at %d kts\n", cd_min, cd_min_kts);
 }
 
-int usage()
+void report(Airplane* a)
 {
-  fprintf(stderr, "Usage: \n");
-  fprintf(stderr, "  yasim <aircraft.xml> [-g [-a meters] [-s kts] [-approach | -cruise] ]\n");
-  fprintf(stderr, "  yasim <aircraft.xml> [-d [-a meters] [-approach | -cruise] ]\n");
-  fprintf(stderr, "  yasim <aircraft.xml> [-m]\n");
-  fprintf(stderr, "                       -g print lift/drag table: aoa, lift, drag, lift/drag \n");
-  fprintf(stderr, "                       -d print drag over TAS: kts, drag\n");
-  fprintf(stderr, "                       -a set altitude in meters!\n");
-  fprintf(stderr, "                       -s set speed in knots\n");
-  fprintf(stderr, "                       -m print mass distribution table: id, x, y, z, mass \n");
-  return 1;
-}
-
-
-int main(int argc, char** argv)
-{
-  FGFDM* fdm = new FGFDM();
-  Airplane* a = fdm->getAirplane();
-
-  if(argc < 2) return usage();
-  // Read
-  try {
-    string file = argv[1];
-    readXML(SGPath(file), *fdm);
-  } 
-  catch (const sg_exception &e) {
-    printf("XML parse error: %s (%s)\n", e.getFormattedMessage().c_str(), e.getOrigin());
-  }
-
-  // ... and run
-  a->compile();
-  if(a->getFailureMsg())
-      printf("SOLUTION FAILURE: %s\n", a->getFailureMsg());
-  if(!a->getFailureMsg() && argc > 2 ) {
-    if(strcmp(argv[2], "-g") == 0) {
-      float alt = 5000, kts = 100;
-      int cfg = CONFIG_NONE;
-      for(int i=3; i<argc; i++) {
-        if (std::strcmp(argv[i], "-a") == 0) {
-          if (i+1 < argc) alt = std::atof(argv[++i]);
-        }
-        else if(std::strcmp(argv[i], "-s") == 0) {
-          if(i+1 < argc) kts = std::atof(argv[++i]);
-        }
-        else if(std::strcmp(argv[i], "-approach") == 0) cfg = CONFIG_APPROACH;
-        else if(std::strcmp(argv[i], "-cruise") == 0) cfg = CONFIG_CRUISE;
-        else return usage();
-      }
-      yasim_graph(a, alt, kts, cfg);
-    } 
-    else if(strcmp(argv[2], "-d") == 0) {
-      float alt = 2000, aoa = a->getCruiseAoA();
-      int cfg = CONFIG_NONE;
-      for(int i=3; i<argc; i++) {
-	if (std::strcmp(argv[i], "-a") == 0) {
-	  if (i+1 < argc) alt = std::atof(argv[++i]);
-	}
-	else if(std::strcmp(argv[i], "-approach") == 0) cfg = CONFIG_APPROACH;
-	else if(std::strcmp(argv[i], "-cruise") == 0) cfg = CONFIG_CRUISE;
-	else return usage();
-      }
-      yasim_drag(a, aoa, alt, cfg);
-    } 
-    else if(strcmp(argv[2], "-m") == 0) {
-      yasim_masses(a);
-    }
-  }
-  else {
     printf("==========================\n");
     printf("= YASim solution results =\n");
     printf("==========================\n");
     float aoa = a->getCruiseAoA() * RAD2DEG;
-    float tailIncidence = -1 * a->getTailIncidence() * RAD2DEG;
+    float tailIncidence = a->getTailIncidence() * RAD2DEG;
     float drag = 1000 * a->getDragCoefficient();
     float cg[3];
     a->getModel()->getBody()->getCG(cg);
@@ -296,6 +229,90 @@ int main(int argc, char** argv)
     printf("  %7.0f, %7.0f, %7.0f\n", SI_inertia[0], SI_inertia[1], SI_inertia[2]);
     printf("  %7.0f, %7.0f, %7.0f\n", SI_inertia[3], SI_inertia[4], SI_inertia[5]);
     printf("  %7.0f, %7.0f, %7.0f\n", SI_inertia[6], SI_inertia[7], SI_inertia[8]);
+}
+
+int usage()
+{
+  fprintf(stderr, "Usage: \n");
+  fprintf(stderr, "  yasim <aircraft.xml> [-g [-a meters] [-s kts] [-approach | -cruise] ]\n");
+  fprintf(stderr, "  yasim <aircraft.xml> [-d [-a meters] [-approach | -cruise] ]\n");
+  fprintf(stderr, "  yasim <aircraft.xml> [-m]\n");
+  fprintf(stderr, "  yasim <aircraft.xml> [-test] [-a meters] [-s kts] [-approach | -cruise] ]\n");
+  fprintf(stderr, "                       -g print lift/drag table: aoa, lift, drag, lift/drag \n");
+  fprintf(stderr, "                       -d print drag over TAS: kts, drag\n");
+  fprintf(stderr, "                       -a set altitude in meters!\n");
+  fprintf(stderr, "                       -s set speed in knots\n");
+  fprintf(stderr, "                       -m print mass distribution table: id, x, y, z, mass \n");
+  fprintf(stderr, "                       -test print summary and output like -g -m \n");
+  return 1;
+}
+
+
+int main(int argc, char** argv)
+{
+  FGFDM* fdm = new FGFDM();
+  Airplane* a = fdm->getAirplane();
+
+  if(argc < 2) return usage();
+  // Read
+  try {
+    string file = argv[1];
+    readXML(SGPath(file), *fdm);
+  } 
+  catch (const sg_exception &e) {
+    printf("XML parse error: %s (%s)\n", e.getFormattedMessage().c_str(), e.getOrigin());
+  }
+
+  // ... and run
+  a->compile();
+  if(a->getFailureMsg())
+      printf("SOLUTION FAILURE: %s\n", a->getFailureMsg());
+  if(!a->getFailureMsg() && argc > 2 ) {
+    bool test = (strcmp(argv[2], "-test") == 0);
+    if((strcmp(argv[2], "-g") == 0) || test) 
+    {
+      float alt = 5000, kts = 100;
+      int cfg = CONFIG_NONE;
+      for(int i=3; i<argc; i++) {
+        if (std::strcmp(argv[i], "-a") == 0) {
+          if (i+1 < argc) alt = std::atof(argv[++i]);
+        }
+        else if(std::strcmp(argv[i], "-s") == 0) {
+          if(i+1 < argc) kts = std::atof(argv[++i]);
+        }
+        else if(std::strcmp(argv[i], "-approach") == 0) cfg = CONFIG_APPROACH;
+        else if(std::strcmp(argv[i], "-cruise") == 0) cfg = CONFIG_CRUISE;
+        else return usage();
+      }
+      if (test) {
+        report(a);
+        printf("\n#-- lift, drag at altitude %.0f meters, %.0f knots, Config %d --\n", alt, kts, cfg);
+      }
+      yasim_graph(a, alt, kts, cfg);
+      if (test) {
+        printf("\n#-- mass distribution --\n");
+        yasim_masses(a);
+      }
+    } 
+    else if(strcmp(argv[2], "-d") == 0) {
+      float alt = 2000, aoa = a->getCruiseAoA();
+      int cfg = CONFIG_NONE;
+      for(int i=3; i<argc; i++) {
+	if (std::strcmp(argv[i], "-a") == 0) {
+	  if (i+1 < argc) alt = std::atof(argv[++i]);
+	}
+	else if(std::strcmp(argv[i], "-approach") == 0) cfg = CONFIG_APPROACH;
+	else if(std::strcmp(argv[i], "-cruise") == 0) cfg = CONFIG_CRUISE;
+	else return usage();
+      }
+      yasim_drag(a, aoa, alt, cfg);
+    } 
+    else if(strcmp(argv[2], "-m") == 0) {
+      yasim_masses(a);
+    }
+  }
+  else {
+    report(a);
   }
   delete fdm;
   return 0;