]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/YASim/Airplane.cpp
YASim now supports the new fuel.nas fuel management system. It
[flightgear.git] / src / FDM / YASim / Airplane.cpp
index 9f1f1de0dbb8f7f12c9583e7b3254e131ba801ea..166a861b54b1b352792da5b002dcfb2ef86d14f2 100644 (file)
@@ -5,11 +5,29 @@
 #include "Glue.hpp"
 #include "RigidBody.hpp"
 #include "Surface.hpp"
+#include "Rotorpart.hpp"
+#include "Rotorblade.hpp"
 #include "Thruster.hpp"
-
 #include "Airplane.hpp"
+
 namespace yasim {
 
+// gadgets
+inline float norm(float f) { return f<1 ? 1/f : f; }
+inline float abs(float f) { return f<0 ? -f : f; }
+
+// Solver threshold.  How close to the solution are we trying
+// to get?  Trying too hard can result in oscillations about
+// the correct solution, which is bad.  Stick this in as a
+// compile time constant for now, and consider making it
+// settable per-model.
+const float STHRESH = 1;
+
+// How slowly do we change values in the solver.  Too slow, and
+// the solution converges very slowly.  Too fast, and it can
+// oscillate.
+const float SOLVE_TWEAK = 0.3226;
+
 Airplane::Airplane()
 {
     _emptyWeight = 0;
@@ -46,13 +64,26 @@ Airplane::~Airplane()
        delete (GearRec*)_gears.get(i);
     for(i=0; i<_surfs.size(); i++)
        delete (Surface*)_surfs.get(i);    
+    for(i=0; i<_contacts.size(); i++)
+        delete[] (float*)_contacts.get(i);
+    for(i=0; i<_solveWeights.size(); i++)
+        delete[] (SolveWeight*)_solveWeights.get(i);
 }
 
 void Airplane::iterate(float dt)
 {
+    // The gear might have moved.  Change their aerodynamics.
+    updateGearState();
+
     _model.iterate();
+}
 
-    // FIXME: Consume fuel
+void Airplane::calcFuelWeights()
+{
+    for(int i=0; i<_tanks.size(); i++) {
+        Tank* t = (Tank*)_tanks.get(i);
+        _model.getBody()->setMass(t->handle, t->fill);
+    }
 }
 
 ControlMap* Airplane::getControlMap()
@@ -71,7 +102,7 @@ void Airplane::getPilotAccel(float* out)
 
     // Gravity
     Glue::geodUp(s->pos, out);
-    Math::mul3(-9.8, out, out);
+    Math::mul3(-9.8f, out, out);
 
     // The regular acceleration
     float tmp[3];
@@ -106,54 +137,42 @@ Gear* Airplane::getGear(int g)
     return ((GearRec*)_gears.get(g))->gear;
 }
 
-void Airplane::setGearState(bool down, float dt)
+void Airplane::updateGearState()
 {
-    int i;
-    for(i=0; i<_gears.size(); i++) {
+    for(int i=0; i<_gears.size(); i++) {
         GearRec* gr = (GearRec*)_gears.get(i);
-        if(gr->time == 0) {
-            // Non-extensible
-            gr->gear->setExtension(1);
-            gr->surf->setXDrag(1);
-            gr->surf->setYDrag(1);
-            gr->surf->setZDrag(1);
-            continue;
-        }
-
-        float diff = dt / gr->time;
-        if(!down) diff = -diff;
-        float ext = gr->gear->getExtension() + diff;
-        if(ext < 0) ext = 0;
-        if(ext > 1) ext = 1;
+        float ext = gr->gear->getExtension();
 
-        gr->gear->setExtension(ext);
         gr->surf->setXDrag(ext);
         gr->surf->setYDrag(ext);
         gr->surf->setZDrag(ext);
     }
 }
 
-void Airplane::setApproach(float speed, float altitude)
-{
-    // The zero AoA will become a calculated stall AoA in compile()
-    setApproach(speed, altitude, 0);
-}
-
-void Airplane::setApproach(float speed, float altitude, float aoa)
+void Airplane::setApproach(float speed, float altitude, float aoa, float fuel)
 {
     _approachSpeed = speed;
     _approachP = Atmosphere::getStdPressure(altitude);
     _approachT = Atmosphere::getStdTemperature(altitude);
     _approachAoA = aoa;
+    _approachFuel = fuel;
 }
  
-void Airplane::setCruise(float speed, float altitude)
+void Airplane::setCruise(float speed, float altitude, float fuel)
 {
     _cruiseSpeed = speed;
     _cruiseP = Atmosphere::getStdPressure(altitude);
     _cruiseT = Atmosphere::getStdTemperature(altitude);
     _cruiseAoA = 0;
     _tailIncidence = 0;
+    _cruiseFuel = fuel;
+}
+
+void Airplane::setElevatorControl(int control)
+{
+    _approachElevator.control = control;
+    _approachElevator.val = 0;
+    _approachControls.add(&_approachElevator);
 }
 
 void Airplane::addApproachControl(int control, float val)
@@ -172,6 +191,15 @@ void Airplane::addCruiseControl(int control, float val)
     _cruiseControls.add(c);
 }
 
+void Airplane::addSolutionWeight(bool approach, int idx, float wgt)
+{
+    SolveWeight* w = new SolveWeight();
+    w->approach = approach;
+    w->idx = idx;
+    w->wgt = wgt;
+    _solveWeights.add(w);
+}
+
 int Airplane::numTanks()
 {
     return _tanks.size();
@@ -182,11 +210,21 @@ float Airplane::getFuel(int tank)
     return ((Tank*)_tanks.get(tank))->fill;
 }
 
+float Airplane::setFuel(int tank, float fuel)
+{
+    ((Tank*)_tanks.get(tank))->fill = fuel;
+}
+
 float Airplane::getFuelDensity(int tank)
 {
     return ((Tank*)_tanks.get(tank))->density;
 }
 
+float Airplane::getTankCapacity(int tank)
+{
+    return ((Tank*)_tanks.get(tank))->cap;
+}
+
 void Airplane::setWeight(float weight)
 {
     _emptyWeight = weight;
@@ -207,7 +245,13 @@ void Airplane::addVStab(Wing* vstab)
     _vstabs.add(vstab);
 }
 
-void Airplane::addFuselage(float* front, float* back, float width)
+void Airplane::addRotor(Rotor* rotor)
+{
+    _rotors.add(rotor);
+}
+
+void Airplane::addFuselage(float* front, float* back, float width,
+                           float taper, float mid)
 {
     Fuselage* f = new Fuselage();
     int i;
@@ -216,6 +260,8 @@ void Airplane::addFuselage(float* front, float* back, float width)
        f->back[i]  = back[i];
     }
     f->width = width;
+    f->taper = taper;
+    f->mid = mid;
     _fuselages.add(f);
 }
 
@@ -231,12 +277,11 @@ int Airplane::addTank(float* pos, float cap, float density)
     return _tanks.add(t);
 }
 
-void Airplane::addGear(Gear* gear, float transitionTime)
+void Airplane::addGear(Gear* gear)
 {
     GearRec* g = new GearRec();
     g->gear = gear;
     g->surf = 0;
-    g->time = transitionTime;
     _gears.add(g);
 }
 
@@ -294,6 +339,7 @@ void Airplane::setFuelFraction(float frac)
     int i;
     for(i=0; i<_tanks.size(); i++) {
         Tank* t = (Tank*)_tanks.get(i);
+        t->fill = frac * t->cap;
         _model.getBody()->setMass(t->handle, t->cap * frac);
     }
 }
@@ -347,8 +393,26 @@ void Airplane::setupState(float aoa, float speed, State* s)
     s->pos[2] = 1;
 }
 
+void Airplane::addContactPoint(float* pos)
+{
+    float* cp = new float[3];
+    cp[0] = pos[0];
+    cp[1] = pos[1];
+    cp[2] = pos[2];
+    _contacts.add(cp);
+}
+
 float Airplane::compileWing(Wing* w)
 {
+    // The tip of the wing is a contact point
+    float tip[3];
+    w->getTip(tip);
+    addContactPoint(tip);
+    if(w->isMirrored()) {
+        tip[1] *= -1;
+        addContactPoint(tip);
+    }
+
     // Make sure it's initialized.  The surfaces will pop out with
     // total drag coefficients equal to their areas, which is what we
     // want.
@@ -358,18 +422,65 @@ float Airplane::compileWing(Wing* w)
     int i;
     for(i=0; i<w->numSurfaces(); i++) {
         Surface* s = (Surface*)w->getSurface(i);
+
+       float td = s->getTotalDrag();
+       s->setTotalDrag(td);
+
         _model.addSurface(s);
 
+        float mass = w->getSurfaceWeight(i);
+        mass = mass * Math::sqrt(mass);
         float pos[3];
         s->getPosition(pos);
-        _model.getBody()->addMass(w->getSurfaceWeight(i), pos);
-        wgt += w->getSurfaceWeight(i);
+        _model.getBody()->addMass(mass, pos);
+        wgt += mass;
+    }
+    return wgt;
+}
+
+float Airplane::compileRotor(Rotor* r)
+{
+    // Todo: add rotor to model!!!
+    // Todo: calc and add mass!!!
+    r->compile();
+    _model.addRotor(r);
+
+    float wgt = 0;
+    int i;
+    for(i=0; i<r->numRotorparts(); i++) {
+        Rotorpart* s = (Rotorpart*)r->getRotorpart(i);
+
+        _model.addRotorpart(s);
+        
+        float mass = s->getWeight();
+        mass = mass * Math::sqrt(mass);
+        float pos[3];
+        s->getPosition(pos);
+        _model.getBody()->addMass(mass, pos);
+        wgt += mass;
+    }
+    
+    for(i=0; i<r->numRotorblades(); i++) {
+        Rotorblade* b = (Rotorblade*)r->getRotorblade(i);
+
+        _model.addRotorblade(b);
+        
+        float mass = b->getWeight();
+        mass = mass * Math::sqrt(mass);
+        float pos[3];
+        b->getPosition(pos);
+        _model.getBody()->addMass(mass, pos);
+        wgt += mass;
     }
     return wgt;
 }
 
 float Airplane::compileFuselage(Fuselage* f)
 {
+    // The front and back are contact points
+    addContactPoint(f->front);
+    addContactPoint(f->back);
+
     float wgt = 0;
     float fwd[3];
     Math::sub3(f->front, f->back, fwd);
@@ -379,12 +490,23 @@ float Airplane::compileFuselage(Fuselage* f)
     float segWgt = len*wid/segs;
     int j;
     for(j=0; j<segs; j++) {
-        float frac = (j+0.5) / segs;
+        float frac = (j+0.5f) / segs;
+
+        float scale = 1;
+        if(frac < f->mid)
+            scale = f->taper+(1-f->taper) * (frac / f->mid);
+        else
+            scale = f->taper+(1-f->taper) * (frac - f->mid) / (1 - f->mid);
+
+        // Where are we?
         float pos[3];
         Math::mul3(frac, fwd, pos);
         Math::add3(f->back, pos, pos);
-        _model.getBody()->addMass(segWgt, pos);
-        wgt += segWgt;
+
+        // _Mass_ weighting goes as surface area^(3/2)
+        float mass = scale*segWgt * Math::sqrt(scale*segWgt);
+        _model.getBody()->addMass(mass, pos);
+        wgt += mass;
 
         // Make a Surface too
         Surface* s = new Surface();
@@ -392,7 +514,7 @@ float Airplane::compileFuselage(Fuselage* f)
        float sideDrag = len/wid;
         s->setYDrag(sideDrag);
         s->setZDrag(sideDrag);
-        s->setTotalDrag(segWgt);
+        s->setTotalDrag(scale*segWgt);
 
         // FIXME: fails for fuselages aligned along the Y axis
         float o[9];
@@ -400,6 +522,8 @@ float Airplane::compileFuselage(Fuselage* f)
         Math::unit3(fwd, x);
         y[0] = 0; y[1] = 1; y[2] = 0;
         Math::cross3(x, y, z);
+       Math::unit3(z, z);
+       Math::cross3(z, x, y);
         s->setOrientation(o);
 
         _model.addSurface(s);
@@ -436,6 +560,39 @@ void Airplane::compileGear(GearRec* gr)
     _surfs.add(s);
 }
 
+void Airplane::compileContactPoints()
+{
+    // Figure it will compress by 20cm
+    float comp[3];
+    float DIST = 0.2f;
+    comp[0] = 0; comp[1] = 0; comp[2] = DIST;
+
+    // Give it a spring constant such that at full compression it will
+    // hold up 10 times the planes mass.  That's about right.  Yeah.
+    float mass = _model.getBody()->getTotalMass();
+    float spring = (1/DIST) * 9.8f * 10.0f * mass;
+    float damp = 2 * Math::sqrt(spring * mass);
+
+    int i;
+    for(i=0; i<_contacts.size(); i++) {
+        float *cp = (float*)_contacts.get(i);
+
+        Gear* g = new Gear();
+        g->setPosition(cp);
+        
+        g->setCompression(comp);
+        g->setSpring(spring);
+        g->setDamping(damp);
+        g->setBrake(1);
+
+        // I made these up
+        g->setStaticFriction(0.6f);
+        g->setDynamicFriction(0.5f);
+
+        _model.addGear(g);
+    }
+}
+
 void Airplane::compile()
 {
     double ground[3];
@@ -451,17 +608,19 @@ void Airplane::compile()
     float aeroWgt = 0;
 
     // The Wing objects
-    aeroWgt += compileWing(_wing);
-    aeroWgt += compileWing(_tail);
+    if (_wing)
+      aeroWgt += compileWing(_wing);
+    if (_tail)
+      aeroWgt += compileWing(_tail);
     int i;
-    for(i=0; i<_vstabs.size(); i++) {
+    for(i=0; i<_vstabs.size(); i++)
         aeroWgt += compileWing((Wing*)_vstabs.get(i)); 
-    }
+    for(i=0; i<_rotors.size(); i++)
+        aeroWgt += compileRotor((Rotor*)_rotors.get(i)); 
     
     // The fuselage(s)
-    for(i=0; i<_fuselages.size(); i++) {
+    for(i=0; i<_fuselages.size(); i++)
         aeroWgt += compileFuselage((Fuselage*)_fuselages.get(i));
-    }
 
     // Count up the absolute weight we have
     float nonAeroWgt = _ballast;
@@ -486,8 +645,8 @@ void Airplane::compile()
         t->handle = body->addMass(0, t->pos);
         totalFuel += t->cap;
     }
-    _cruiseWeight = _emptyWeight + totalFuel*0.5;
-    _approachWeight = _emptyWeight + totalFuel*0.2;
+    _cruiseWeight = _emptyWeight + totalFuel*0.5f;
+    _approachWeight = _emptyWeight + totalFuel*0.2f;
 
     body->recalc();
 
@@ -503,14 +662,18 @@ void Airplane::compile()
 
     // Ground effect
     float gepos[3];
-    float gespan = _wing->getGroundEffect(gepos);
-    _model.setGroundEffect(gepos, gespan, .3);
+    float gespan = 0;
+    if(_wing)
+      gespan = _wing->getGroundEffect(gepos);
+    _model.setGroundEffect(gepos, gespan, 0.15f);
 
     solveGear();
-    solve();
+    if(_wing && _tail) solve();
+    else solveHelicopter();
 
-    // Drop the gear (use a really big dt)
-    setGearState(true, 1000000);
+    // Do this after solveGear, because it creates "gear" objects that
+    // we don't want to affect.
+    compileContactPoints();
 }
 
 void Airplane::solveGear()
@@ -530,7 +693,7 @@ void Airplane::solveGear()
         Gear* g = gr->gear;
         g->getPosition(pos);
        Math::sub3(cg, pos, pos);
-        gr->wgt = 1/(0.5+Math::sqrt(pos[0]*pos[0] + pos[1]*pos[1]));
+        gr->wgt = 1.0f/(0.5f+Math::sqrt(pos[0]*pos[0] + pos[1]*pos[1]));
         total += gr->wgt;
     }
 
@@ -539,14 +702,14 @@ void Airplane::solveGear()
         ((GearRec*)_gears.get(i))->wgt /= total;
     
     // The force at max compression should be sufficient to stop a
-    // plane moving downwards at 3x the approach descent rate.  Assume
+    // plane moving downwards at 2x the approach descent rate.  Assume
     // a 3 degree approach.
-    float descentRate = 3*_approachSpeed/19.1;
+    float descentRate = 2.0f*_approachSpeed/19.1f;
 
     // Spread the kinetic energy according to the gear weights.  This
     // will results in an equal compression fraction (not distance) of
     // each gear.
-    float energy = 0.5*_approachWeight*descentRate*descentRate;
+    float energy = 0.5f*_approachWeight*descentRate*descentRate;
 
     for(i=0; i<_gears.size(); i++) {
         GearRec* gr = (GearRec*)_gears.get(i);
@@ -558,14 +721,23 @@ void Airplane::solveGear()
         // Energy in a spring: e = 0.5 * k * len^2
         float k = 2 * e / (len*len);
 
-        gr->gear->setSpring(k);
+        gr->gear->setSpring(k * gr->gear->getSpring());
 
         // Critically damped (too damped, too!)
-        gr->gear->setDamping(2*Math::sqrt(k*_approachWeight*gr->wgt));
+        gr->gear->setDamping(2*Math::sqrt(k*_approachWeight*gr->wgt)
+                             * gr->gear->getDamping());
 
         // These are pretty generic
-        gr->gear->setStaticFriction(0.8);
-        gr->gear->setDynamicFriction(0.7);
+        gr->gear->setStaticFriction(0.8f);
+        gr->gear->setDynamicFriction(0.7f);
+    }
+}
+
+void Airplane::initEngines()
+{
+    for(int i=0; i<_thrusters.size(); i++) {
+        ThrustRec* tr = (ThrustRec*)_thrusters.get(i);
+       tr->thruster->init();
     }
 }
 
@@ -576,11 +748,24 @@ void Airplane::stabilizeThrust()
        _model.getThruster(i)->stabilize();
 }
 
+void Airplane::setupWeights(bool isApproach)
+{
+    int i;
+    for(i=0; i<_weights.size(); i++)
+        setWeight(i, 0);
+    for(i=0; i<_solveWeights.size(); i++) {
+        SolveWeight* w = (SolveWeight*)_solveWeights.get(i);
+        if(w->approach == isApproach)
+            setWeight(w->idx, w->wgt);
+    }
+}
+
 void Airplane::runCruise()
 {
     setupState(_cruiseAoA, _cruiseSpeed, &_cruiseState);
     _model.setState(&_cruiseState);
-    _model.setAir(_cruiseP, _cruiseT);
+    _model.setAir(_cruiseP, _cruiseT,
+                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
 
     // The control configuration
     _controls.reset();
@@ -589,29 +774,30 @@ void Airplane::runCruise()
        Control* c = (Control*)_cruiseControls.get(i);
        _controls.setInput(c->control, c->val);
     }
-    _controls.applyControls();
+    _controls.applyControls(1000000); // Huge dt value
 
     // The local wind
     float wind[3];
     Math::mul3(-1, _cruiseState.v, wind);
     Math::vmul33(_cruiseState.orient, wind, wind);
  
-    // Gear are up (if they're non-retractable, this is a noop)
-    setGearState(false, 100000);
-    
-    // Cruise is by convention at 50% tank capacity
-    setFuelFraction(0.5);
+    setFuelFraction(_cruiseFuel);
+    setupWeights(false);
    
     // Set up the thruster parameters and iterate until the thrust
     // stabilizes.
     for(i=0; i<_thrusters.size(); i++) {
        Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
        t->setWind(wind);
-       t->setAir(_cruiseP, _cruiseT);
+       t->setAir(_cruiseP, _cruiseT,
+                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
     }
     stabilizeThrust();
 
+    updateGearState();
+
     // Precompute thrust in the model, and calculate aerodynamic forces
+    _model.getBody()->recalc();
     _model.getBody()->reset();
     _model.initIteration();
     _model.calcForces(&_cruiseState);
@@ -621,7 +807,8 @@ void Airplane::runApproach()
 {
     setupState(_approachAoA, _approachSpeed, &_approachState);
     _model.setState(&_approachState);
-    _model.setAir(_approachP, _approachT);
+    _model.setAir(_approachP, _approachT,
+                  Atmosphere::calcStdDensity(_approachP, _approachT));
 
     // The control configuration
     _controls.reset();
@@ -630,29 +817,31 @@ void Airplane::runApproach()
        Control* c = (Control*)_approachControls.get(i);
        _controls.setInput(c->control, c->val);
     }
-    _controls.applyControls();
+    _controls.applyControls(1000000);
 
     // The local wind
     float wind[3];
     Math::mul3(-1, _approachState.v, wind);
     Math::vmul33(_approachState.orient, wind, wind);
     
-    // Approach is by convention at 20% tank capacity
-    setFuelFraction(0.2);
+    setFuelFraction(_approachFuel);
 
-    // Gear are down
-    setGearState(true, 100000);
+    setupWeights(true);
 
     // Run the thrusters until they get to a stable setting.  FIXME:
     // this is lots of wasted work.
     for(i=0; i<_thrusters.size(); i++) {
        Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
        t->setWind(wind);
-       t->setAir(_approachP, _approachT);
+       t->setAir(_approachP, _approachT,
+                  Atmosphere::calcStdDensity(_approachP, _approachT));
     }
     stabilizeThrust();
 
+    updateGearState();
+
     // Precompute thrust in the model, and calculate aerodynamic forces
+    _model.getBody()->recalc();
     _model.getBody()->reset();
     _model.initIteration();
     _model.calcForces(&_approachState);
@@ -660,10 +849,12 @@ void Airplane::runApproach()
 
 void Airplane::applyDragFactor(float factor)
 {
-    float applied = Math::sqrt(factor);
+    float applied = Math::pow(factor, SOLVE_TWEAK);
     _dragFactor *= applied;
-    _wing->setDragScale(_wing->getDragScale() * applied);
-    _tail->setDragScale(_tail->getDragScale() * applied);
+    if(_wing)
+      _wing->setDragScale(_wing->getDragScale() * applied);
+    if(_tail)
+      _tail->setDragScale(_tail->getDragScale() * applied);
     int i;
     for(i=0; i<_vstabs.size(); i++) {
        Wing* w = (Wing*)_vstabs.get(i);
@@ -677,10 +868,12 @@ void Airplane::applyDragFactor(float factor)
 
 void Airplane::applyLiftRatio(float factor)
 {
-    float applied = Math::sqrt(factor);
+    float applied = Math::pow(factor, SOLVE_TWEAK);
     _liftRatio *= applied;
-    _wing->setLiftRatio(_wing->getLiftRatio() * applied);
-    _tail->setLiftRatio(_tail->getLiftRatio() * applied);
+    if(_wing)
+      _wing->setLiftRatio(_wing->getLiftRatio() * applied);
+    if(_tail)
+      _tail->setLiftRatio(_tail->getLiftRatio() * applied);
     int i;
     for(i=0; i<_vstabs.size(); i++) {
         Wing* w = (Wing*)_vstabs.get(i);
@@ -704,15 +897,16 @@ float Airplane::normFactor(float f)
 
 void Airplane::solve()
 {
-    static const float ARCMIN = 0.0002909;
+    static const float ARCMIN = 0.0002909f;
 
     float tmp[3];
     _solutionIterations = 0;
     _failureMsg = 0;
+
     while(1) {
-       if(_solutionIterations++ > 10000) {
+        if(_solutionIterations++ > 10000) { 
             _failureMsg = "Solution failed to converge after 10000 iterations";
-           return;
+            return;
         }
 
        // Run an iteration at cruise, and extract the needed numbers:
@@ -722,16 +916,23 @@ void Airplane::solve()
        float thrust = tmp[0];
 
        _model.getBody()->getAccel(tmp);
+        Math::tmul33(_cruiseState.orient, tmp, tmp);
        float xforce = _cruiseWeight * tmp[0];
        float clift0 = _cruiseWeight * tmp[2];
 
        _model.getBody()->getAngularAccel(tmp);
+        Math::tmul33(_cruiseState.orient, tmp, tmp);
        float pitch0 = tmp[1];
 
        // Run an approach iteration, and do likewise
        runApproach();
 
+       _model.getBody()->getAngularAccel(tmp);
+        Math::tmul33(_approachState.orient, tmp, tmp);
+       double apitch0 = tmp[1];
+
        _model.getBody()->getAccel(tmp);
+        Math::tmul33(_approachState.orient, tmp, tmp);
        float alift = _approachWeight * tmp[2];
 
        // Modify the cruise AoA a bit to get a derivative
@@ -740,6 +941,7 @@ void Airplane::solve()
        _cruiseAoA -= ARCMIN;
 
        _model.getBody()->getAccel(tmp);
+        Math::tmul33(_cruiseState.orient, tmp, tmp);
        float clift1 = _cruiseWeight * tmp[2];
 
        // Do the same with the tail incidence
@@ -748,10 +950,11 @@ void Airplane::solve()
        _tail->setIncidence(_tailIncidence);
 
        _model.getBody()->getAngularAccel(tmp);
+        Math::tmul33(_cruiseState.orient, tmp, tmp);
        float pitch1 = tmp[1];
 
        // Now calculate:
-       float awgt = 9.8 * _approachWeight;
+       float awgt = 9.8f * _approachWeight;
 
        float dragFactor = thrust / (thrust-xforce);
        float liftFactor = awgt / (awgt+alift);
@@ -759,36 +962,59 @@ void Airplane::solve()
        float tailDelta = -pitch0 * (ARCMIN/(pitch1-pitch0));
 
         // Sanity:
-        if(dragFactor <= 0) {
-            _failureMsg = "Zero or negative drag adjustment.";
-            return;
-        } else if(liftFactor <= 0) {
-            _failureMsg = "Zero or negative lift adjustment.";
-            return;
-        }
+        if(dragFactor <= 0 || liftFactor <= 0)
+            break;
+
+        // And the elevator control in the approach.  This works just
+        // like the tail incidence computation (it's solving for the
+        // same thing -- pitching moment -- by diddling a different
+        // variable).
+        const float ELEVDIDDLE = 0.001f;
+        _approachElevator.val += ELEVDIDDLE;
+        runApproach();
+        _approachElevator.val -= ELEVDIDDLE;
+
+       _model.getBody()->getAngularAccel(tmp);
+        Math::tmul33(_approachState.orient, tmp, tmp);
+       double apitch1 = tmp[1];
+        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
+        // numbers in the right ballpark.
 
-       // And apply:
        applyDragFactor(dragFactor);
        applyLiftRatio(liftFactor);
 
        // DON'T do the following until the above are sane
-       if(normFactor(dragFactor) > 1.1
-          || normFactor(liftFactor) > 1.1)
+       if(normFactor(dragFactor) > STHRESH*1.0001
+          || normFactor(liftFactor) > STHRESH*1.0001)
        {
            continue;
        }
 
-       // OK, now we can adjust the minor variables
-       _cruiseAoA += 0.5*aoaDelta;
-       _tailIncidence += 0.5*tailDelta;
+       // OK, now we can adjust the minor variables:
+       _cruiseAoA += SOLVE_TWEAK*aoaDelta;
+       _tailIncidence += SOLVE_TWEAK*tailDelta;
        
-       _cruiseAoA = clamp(_cruiseAoA, -.174, .174);
-       _tailIncidence = clamp(_tailIncidence, -.174, .174);
+       _cruiseAoA = clamp(_cruiseAoA, -0.175f, 0.175f);
+       _tailIncidence = clamp(_tailIncidence, -0.175f, 0.175f);
 
-        if(dragFactor < 1.00001 && liftFactor < 1.00001 &&
-           aoaDelta < .000017   && tailDelta < .000017)
+        if(abs(xforce/_cruiseWeight) < STHRESH*0.0001 &&
+           abs(alift/_approachWeight) < STHRESH*0.0001 &&
+           abs(aoaDelta) < STHRESH*.000017 &&
+           abs(tailDelta) < STHRESH*.000017)
         {
-            break;
+            // 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";
+                break;
+            }
         }
     }
 
@@ -798,12 +1024,26 @@ void Airplane::solve()
     } else if(_liftRatio < 1e-04 || _liftRatio > 1e4) {
        _failureMsg = "Lift ratio beyond reasonable bounds.";
        return;
-    } else if(Math::abs(_cruiseAoA) >= .174) {
+    } else if(Math::abs(_cruiseAoA) >= .17453293) {
        _failureMsg = "Cruise AoA > 10 degrees";
        return;
-    } else if(Math::abs(_tailIncidence) >= .174) {
+    } else if(Math::abs(_tailIncidence) >= .17453293) {
        _failureMsg = "Tail incidence > 10 degrees";
        return;
     }
 }
+
+void Airplane::solveHelicopter()
+{
+    _solutionIterations = 0;
+    _failureMsg = 0;
+
+    applyDragFactor(Math::pow(15.7/1000, 1/SOLVE_TWEAK));
+    applyLiftRatio(Math::pow(104, 1/SOLVE_TWEAK));
+    setupState(0,0, &_cruiseState);
+    _model.setState(&_cruiseState);
+    _controls.reset();
+    _model.getBody()->reset();
+}
+
 }; // namespace yasim