]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/YASim/Airplane.cpp
FGPUIDialog: fix reading from already free'd memory.
[flightgear.git] / src / FDM / YASim / Airplane.cpp
index 7132772aaad2c0b3a974fd18538396831c312107..15fd4d683bb1d1acf3096425bf379ff9a59f4a6e 100644 (file)
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#  include "config.h"
+#endif
+
 #include "Atmosphere.hpp"
 #include "ControlMap.hpp"
 #include "Gear.hpp"
@@ -5,11 +9,29 @@
 #include "Glue.hpp"
 #include "RigidBody.hpp"
 #include "Surface.hpp"
+#include "Rotorpart.hpp"
 #include "Thruster.hpp"
-
+#include "Hitch.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;
@@ -21,16 +43,20 @@ Airplane::Airplane()
     _cruiseT = 0;
     _cruiseSpeed = 0;
     _cruiseWeight = 0;
+    _cruiseGlideAngle = 0;
     _approachP = 0;
     _approachT = 0;
     _approachSpeed = 0;
     _approachAoA = 0;
     _approachWeight = 0;
+    _approachGlideAngle = 0;
 
     _dragFactor = 1;
     _liftRatio = 1;
     _cruiseAoA = 0;
     _tailIncidence = 0;
+
+    _failureMsg = 0;
 }
 
 Airplane::~Airplane()
@@ -42,19 +68,49 @@ Airplane::~Airplane()
        delete (Tank*)_tanks.get(i);
     for(i=0; i<_thrusters.size(); i++)
        delete (ThrustRec*)_thrusters.get(i);
-    for(i=0; i<_gears.size(); i++)
-       delete (GearRec*)_gears.get(i);
+    for(i=0; i<_gears.size(); i++) {
+       GearRec* g = (GearRec*)_gears.get(i);
+        delete g->gear;
+        delete g;
+    }
     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<_contacts.size(); i++) {
+        ContactRec* c = (ContactRec*)_contacts.get(i);
+        delete c->gear;
+        delete c;
+    }
+    for(i=0; i<_solveWeights.size(); i++)
+        delete (SolveWeight*)_solveWeights.get(i);
+    for(i=0; i<_cruiseControls.size(); i++)
+        delete (Control*)_cruiseControls.get(i);
+    for(i=0; i<_approachControls.size(); i++) {
+        Control* c = (Control*)_approachControls.get(i);
+        if(c != &_approachElevator)
+            delete c;
+    }
+    delete _wing;
+    delete _tail;
+    for(i=0; i<_vstabs.size(); i++)
+        delete (Wing*)_vstabs.get(i);
+    for(i=0; i<_weights.size(); i++)
+        delete (WeightRec*)_weights.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()
@@ -73,15 +129,18 @@ void Airplane::getPilotAccel(float* out)
 
     // Gravity
     Glue::geodUp(s->pos, out);
-    Math::mul3(-9.8, out, out);
+    Math::mul3(-9.8f, out, out);
+    Math::vmul33(s->orient, out, out);
+    out[0] = -out[0];
 
     // The regular acceleration
     float tmp[3];
-    Math::mul3(-1, s->acc, tmp);
-    Math::add3(tmp, out, out);
-
     // Convert to aircraft coordinates
-    Math::vmul33(s->orient, out, out);
+    Math::vmul33(s->orient, s->acc, tmp);
+    tmp[1] = -tmp[1];
+    tmp[2] = -tmp[2];
+
+    Math::add3(tmp, out, out);
 
     // FIXME: rotational & centripetal acceleration needed
 }
@@ -108,54 +167,59 @@ Gear* Airplane::getGear(int g)
     return ((GearRec*)_gears.get(g))->gear;
 }
 
-void Airplane::setGearState(bool down, float dt)
+Hook* Airplane::getHook()
 {
-    int i;
-    for(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;
-        }
+    return _model.getHook();
+}
+
+Launchbar* Airplane::getLaunchbar()
+{
+    return _model.getLaunchbar();
+}
+
+Rotorgear* Airplane::getRotorgear()
+{
+    return _model.getRotorgear();
+}
 
-        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;
+void Airplane::updateGearState()
+{
+    for(int i=0; i<_gears.size(); i++) {
+        GearRec* gr = (GearRec*)_gears.get(i);
+        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, float gla)
 {
     _approachSpeed = speed;
     _approachP = Atmosphere::getStdPressure(altitude);
     _approachT = Atmosphere::getStdTemperature(altitude);
     _approachAoA = aoa;
+    _approachFuel = fuel;
+    _approachGlideAngle = gla;
 }
  
-void Airplane::setCruise(float speed, float altitude)
+void Airplane::setCruise(float speed, float altitude, float fuel, float gla)
 {
     _cruiseSpeed = speed;
     _cruiseP = Atmosphere::getStdPressure(altitude);
     _cruiseT = Atmosphere::getStdTemperature(altitude);
     _cruiseAoA = 0;
     _tailIncidence = 0;
+    _cruiseFuel = fuel;
+    _cruiseGlideAngle = gla;
+}
+
+void Airplane::setElevatorControl(int control)
+{
+    _approachElevator.control = control;
+    _approachElevator.val = 0;
+    _approachControls.add(&_approachElevator);
 }
 
 void Airplane::addApproachControl(int control, float val)
@@ -174,6 +238,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();
@@ -184,11 +257,21 @@ float Airplane::getFuel(int tank)
     return ((Tank*)_tanks.get(tank))->fill;
 }
 
+float Airplane::setFuel(int tank, float fuel)
+{
+    return ((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;
@@ -210,7 +293,8 @@ void Airplane::addVStab(Wing* vstab)
 }
 
 void Airplane::addFuselage(float* front, float* back, float width,
-                           float taper, float mid)
+                           float taper, float mid, 
+                           float cx, float cy, float cz, float idrag)
 {
     Fuselage* f = new Fuselage();
     int i;
@@ -221,6 +305,10 @@ void Airplane::addFuselage(float* front, float* back, float width,
     f->width = width;
     f->taper = taper;
     f->mid = mid;
+    f->_cx=cx;
+    f->_cy=cy;
+    f->_cz=cz;
+    f->_idrag=idrag;
     _fuselages.add(f);
 }
 
@@ -236,15 +324,29 @@ 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);
 }
 
+void Airplane::addHook(Hook* hook)
+{
+    _model.addHook(hook);
+}
+
+void Airplane::addHitch(Hitch* hitch)
+{
+    _model.addHitch(hitch);
+}
+
+void Airplane::addLaunchbar(Launchbar* launchbar)
+{
+    _model.addLaunchbar(launchbar);
+}
+
 void Airplane::addThruster(Thruster* thruster, float mass, float* cg)
 {
     ThrustRec* t = new ThrustRec();
@@ -299,6 +401,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);
     }
 }
@@ -323,7 +426,7 @@ float Airplane::getTailIncidence()
     return _tailIncidence;
 }
 
-char* Airplane::getFailureMsg()
+const char* Airplane::getFailureMsg()
 {
     return _failureMsg;
 }
@@ -333,7 +436,7 @@ int Airplane::getSolutionIterations()
     return _solutionIterations;
 }
 
-void Airplane::setupState(float aoa, float speed, State* s)
+void Airplane::setupState(float aoa, float speed, float gla, State* s)
 {
     float cosAoA = Math::cos(aoa);
     float sinAoA = Math::sin(aoa);
@@ -341,7 +444,7 @@ void Airplane::setupState(float aoa, float speed, State* s)
     s->orient[3] =       0; s->orient[4] = 1; s->orient[5] =      0;
     s->orient[6] = -sinAoA; s->orient[7] = 0; s->orient[8] = cosAoA;
 
-    s->v[0] = speed; s->v[1] = 0; s->v[2] = 0;
+    s->v[0] = speed*Math::cos(gla); s->v[1] = -speed*Math::sin(gla); s->v[2] = 0;
 
     int i;
     for(i=0; i<3; i++)
@@ -354,11 +457,12 @@ void Airplane::setupState(float aoa, float speed, State* s)
 
 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);
+    ContactRec* c = new ContactRec;
+    c->gear = 0;
+    c->p[0] = pos[0];
+    c->p[1] = pos[1];
+    c->p[2] = pos[2];
+    _contacts.add(c);
 }
 
 float Airplane::compileWing(Wing* w)
@@ -397,6 +501,11 @@ float Airplane::compileWing(Wing* w)
     return wgt;
 }
 
+void Airplane::compileRotorgear()
+{
+    getRotorgear()->compile();
+}
+
 float Airplane::compileFuselage(Fuselage* f)
 {
     // The front and back are contact points
@@ -407,12 +516,16 @@ float Airplane::compileFuselage(Fuselage* f)
     float fwd[3];
     Math::sub3(f->front, f->back, fwd);
     float len = Math::mag3(fwd);
+    if (len == 0) {
+        _failureMsg = "Zero length fuselage";
+       return 0;
+    }
     float wid = f->width;
     int segs = (int)Math::ceil(len/wid);
     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)
@@ -434,9 +547,10 @@ float Airplane::compileFuselage(Fuselage* f)
         Surface* s = new Surface();
         s->setPosition(pos);
        float sideDrag = len/wid;
-        s->setYDrag(sideDrag);
-        s->setZDrag(sideDrag);
-        s->setTotalDrag(scale*segWgt);
+        s->setYDrag(sideDrag*f->_cy);
+        s->setZDrag(sideDrag*f->_cz);
+        s->setTotalDrag(scale*segWgt*f->_cx);
+        s->setInducedDrag(f->_idrag);
 
         // FIXME: fails for fuselages aligned along the Y axis
         float o[9];
@@ -444,6 +558,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);
@@ -484,21 +600,22 @@ void Airplane::compileContactPoints()
 {
     // Figure it will compress by 20cm
     float comp[3];
-    float DIST = 0.2;
+    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.8 * 10 * mass;
+    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);
+        ContactRec* c = (ContactRec*)_contacts.get(i);
 
         Gear* g = new Gear();
-        g->setPosition(cp);
+        c->gear = g;
+        g->setPosition(c->p);
         
         g->setCompression(comp);
         g->setSpring(spring);
@@ -506,8 +623,9 @@ void Airplane::compileContactPoints()
         g->setBrake(1);
 
         // I made these up
-        g->setStaticFriction(0.6);
-        g->setDynamicFriction(0.5);
+        g->setStaticFriction(0.6f);
+        g->setDynamicFriction(0.5f);
+        g->setContactPoint(1);
 
         _model.addGear(g);
     }
@@ -515,10 +633,6 @@ void Airplane::compileContactPoints()
 
 void Airplane::compile()
 {
-    double ground[3];
-    ground[0] = 0; ground[1] = 0; ground[2] = 1;
-    _model.setGroundPlane(ground, -100000);
-
     RigidBody* body = _model.getBody();
     int firstMass = body->numMasses();
 
@@ -528,17 +642,18 @@ 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)); 
-    }
-    
+
+
     // 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;
@@ -563,8 +678,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*_cruiseFuel;
+    _approachWeight = _emptyWeight + totalFuel*_approachFuel;
 
     body->recalc();
 
@@ -579,19 +694,29 @@ void Airplane::compile()
     }
 
     // Ground effect
-    float gepos[3];
-    float gespan = _wing->getGroundEffect(gepos);
-    _model.setGroundEffect(gepos, gespan, .3);
+    if(_wing) {
+        float gepos[3];
+        float gespan = 0;
+        gespan = _wing->getGroundEffect(gepos);
+        _model.setGroundEffect(gepos, gespan, 0.15f);
+    }
+
+    // solve function below resets failure message
+    // so check if we have any problems and abort here
+    if (_failureMsg) return;
 
     solveGear();
-    solve();
+    if(_wing && _tail) solve();
+    else
+    {
+       // The rotor(s) mass:
+       compileRotorgear(); 
+       solveHelicopter();
+    }
 
     // Do this after solveGear, because it creates "gear" objects that
     // we don't want to affect.
     compileContactPoints();
-
-    // Drop the gear (use a really big dt)
-    setGearState(true, 1000000);
 }
 
 void Airplane::solveGear()
@@ -611,8 +736,9 @@ 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]));
-        total += gr->wgt;
+        gr->wgt = 1.0f/(0.5f+Math::sqrt(pos[0]*pos[0] + pos[1]*pos[1]));
+        if (!g->getIgnoreWhileSolving())
+            total += gr->wgt;
     }
 
     // Renormalize so they sum to 1
@@ -620,33 +746,30 @@ 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);
         float e = energy * gr->wgt;
         float comp[3];
         gr->gear->getCompression(comp);
-        float len = Math::mag3(comp);
+        float len = Math::mag3(comp)*(1+2*gr->gear->getInitialLoad());
 
         // 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));
-
-        // These are pretty generic
-        gr->gear->setStaticFriction(0.8);
-        gr->gear->setDynamicFriction(0.7);
+        gr->gear->setDamping(2*Math::sqrt(k*_approachWeight*gr->wgt)
+                             * gr->gear->getDamping());
     }
 }
 
@@ -665,11 +788,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);
+    setupState(_cruiseAoA, _cruiseSpeed,_cruiseGlideAngle, &_cruiseState);
     _model.setState(&_cruiseState);
-    _model.setAir(_cruiseP, _cruiseT);
+    _model.setAir(_cruiseP, _cruiseT,
+                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
 
     // The control configuration
     _controls.reset();
@@ -678,29 +814,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);
@@ -708,9 +845,10 @@ void Airplane::runCruise()
 
 void Airplane::runApproach()
 {
-    setupState(_approachAoA, _approachSpeed, &_approachState);
+    setupState(_approachAoA, _approachSpeed,_approachGlideAngle, &_approachState);
     _model.setState(&_approachState);
-    _model.setAir(_approachP, _approachT);
+    _model.setAir(_approachP, _approachT,
+                  Atmosphere::calcStdDensity(_approachP, _approachT));
 
     // The control configuration
     _controls.reset();
@@ -719,29 +857,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);
@@ -749,10 +889,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);
@@ -766,10 +908,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);
@@ -793,34 +937,42 @@ 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:
        runCruise();
 
        _model.getThrust(tmp);
-       float thrust = tmp[0];
+        float thrust = tmp[0] + _cruiseWeight * Math::sin(_cruiseGlideAngle) * 9.81;
 
        _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
@@ -829,6 +981,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
@@ -837,10 +990,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);
@@ -848,36 +1002,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;
+            }
         }
     }
 
@@ -887,12 +1064,42 @@ 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;
+    if (getRotorgear()!=0)
+    {
+        Rotorgear* rg = getRotorgear();
+        applyDragFactor(Math::pow(rg->getYasimDragFactor()/1000,
+            1/SOLVE_TWEAK));
+        applyLiftRatio(Math::pow(rg->getYasimLiftFactor(),
+            1/SOLVE_TWEAK));
+    }
+    else
+    //huh, no wing and no rotor? (_rotorgear is constructed, 
+    //if a rotor is defined
+    {
+        applyDragFactor(Math::pow(15.7/1000, 1/SOLVE_TWEAK));
+        applyLiftRatio(Math::pow(104, 1/SOLVE_TWEAK));
+    }
+    setupState(0,0,0, &_cruiseState);
+    _model.setState(&_cruiseState);
+    setupWeights(true);
+    _controls.reset();
+    _model.getBody()->reset();
+    _model.setAir(_cruiseP, _cruiseT,
+                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
+    
+}
+
 }; // namespace yasim