]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/YASim/Airplane.cpp
Oops. Stray checkin. Revert.
[flightgear.git] / src / FDM / YASim / Airplane.cpp
index 3d76fff548560e8259319cf1efaa38c6aa1d18fa..82b10d7d144844c255759b23206c5a6867243f0e 100644 (file)
 #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; }
+
 Airplane::Airplane()
 {
     _emptyWeight = 0;
@@ -35,20 +39,26 @@ Airplane::Airplane()
 
 Airplane::~Airplane()
 {
-    for(int i=0; i<_fuselages.size(); i++)
+    int i;
+    for(i=0; i<_fuselages.size(); i++)
        delete (Fuselage*)_fuselages.get(i);
-    for(int i=0; i<_tanks.size(); i++)
+    for(i=0; i<_tanks.size(); i++)
        delete (Tank*)_tanks.get(i);
-    for(int i=0; i<_thrusters.size(); i++)
+    for(i=0; i<_thrusters.size(); i++)
        delete (ThrustRec*)_thrusters.get(i);
-    for(int i=0; i<_gears.size(); i++)
+    for(i=0; i<_gears.size(); i++)
        delete (GearRec*)_gears.get(i);
-    for(int i=0; i<_surfs.size(); 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);
 }
 
 void Airplane::iterate(float dt)
 {
+    // The gear might have moved.  Change their aerodynamics.
+    updateGearState();
+
     _model.iterate();
 
     // FIXME: Consume fuel
@@ -70,7 +80,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];
@@ -85,12 +95,14 @@ void Airplane::getPilotAccel(float* out)
 
 void Airplane::setPilotPos(float* pos)
 {
-    for(int i=0; i<3; i++) _pilotPos[i] = pos[i];
+    int i;
+    for(i=0; i<3; i++) _pilotPos[i] = pos[i];
 }
 
 void Airplane::getPilotPos(float* out)
 {
-    for(int i=0; i<3; i++) out[i] = _pilotPos[i];
+    int i;
+    for(i=0; i<3; i++) out[i] = _pilotPos[i];
 }
 
 int Airplane::numGear()
@@ -103,26 +115,12 @@ Gear* Airplane::getGear(int g)
     return ((GearRec*)_gears.get(g))->gear;
 }
 
-void Airplane::setGearState(bool down, float dt)
+void Airplane::updateGearState()
 {
     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);
@@ -152,6 +150,13 @@ void Airplane::setCruise(float speed, float altitude)
     _tailIncidence = 0;
 }
 
+void Airplane::setElevatorControl(int control)
+{
+    _approachElevator.control = control;
+    _approachElevator.val = 0;
+    _approachControls.add(&_approachElevator);
+}
+
 void Airplane::addApproachControl(int control, float val)
 {
     Control* c = new Control();
@@ -203,21 +208,26 @@ void Airplane::addVStab(Wing* vstab)
     _vstabs.add(vstab);
 }
 
-void Airplane::addFuselage(float* front, float* back, float width)
+void Airplane::addFuselage(float* front, float* back, float width,
+                           float taper, float mid)
 {
     Fuselage* f = new Fuselage();
-    for(int i=0; i<3; i++) {
+    int i;
+    for(i=0; i<3; i++) {
        f->front[i] = front[i];
        f->back[i]  = back[i];
     }
     f->width = width;
+    f->taper = taper;
+    f->mid = mid;
     _fuselages.add(f);
 }
 
 int Airplane::addTank(float* pos, float cap, float density)
 {
     Tank* t = new Tank();
-    for(int i=0; i<3; i++) t->pos[i] = pos[i];
+    int i;
+    for(i=0; i<3; i++) t->pos[i] = pos[i];
     t->cap = cap;
     t->fill = cap;
     t->density = density;
@@ -225,12 +235,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);
 }
 
@@ -239,7 +248,8 @@ void Airplane::addThruster(Thruster* thruster, float mass, float* cg)
     ThrustRec* t = new ThrustRec();
     t->thruster = thruster;
     t->mass = mass;
-    for(int i=0; i<3; i++) t->cg[i] = cg[i];
+    int i;
+    for(i=0; i<3; i++) t->cg[i] = cg[i];
     _thrusters.add(t);
 }
 
@@ -284,7 +294,8 @@ void Airplane::setWeight(int handle, float mass)
 
 void Airplane::setFuelFraction(float frac)
 {
-    for(int i=0; i<_tanks.size(); i++) {
+    int i;
+    for(i=0; i<_tanks.size(); i++) {
         Tank* t = (Tank*)_tanks.get(i);
         _model.getBody()->setMass(t->handle, t->cap * frac);
     }
@@ -330,7 +341,8 @@ void Airplane::setupState(float aoa, float speed, State* s)
 
     s->v[0] = speed; s->v[1] = 0; s->v[2] = 0;
 
-    for(int i=0; i<3; i++)
+    int i;
+    for(i=0; i<3; i++)
        s->pos[i] = s->rot[i] = s->acc[i] = s->racc[i] = 0;
 
     // Put us 1m above the origin, or else the gravity computation in
@@ -338,28 +350,57 @@ 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.
     w->compile();
 
     float wgt = 0;
-    for(int i=0; i<w->numSurfaces(); i++) {
+    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::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);
@@ -367,13 +408,25 @@ float Airplane::compileFuselage(Fuselage* f)
     float wid = f->width;
     int segs = (int)Math::ceil(len/wid);
     float segWgt = len*wid/segs;
-    for(int j=0; j<segs; j++) {
-        float frac = (j+0.5) / segs;
+    int j;
+    for(j=0; j<segs; j++) {
+        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();
@@ -381,7 +434,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];
@@ -389,6 +442,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);
@@ -425,6 +480,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];
@@ -442,49 +530,50 @@ void Airplane::compile()
     // The Wing objects
     aeroWgt += compileWing(_wing);
     aeroWgt += compileWing(_tail);
-    for(int i=0; i<_vstabs.size(); i++) {
+    int i;
+    for(i=0; i<_vstabs.size(); i++) {
         aeroWgt += compileWing((Wing*)_vstabs.get(i)); 
     }
     
     // The fuselage(s)
-    for(int 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;
-    for(int i=0; i<_thrusters.size(); i++)
+    for(i=0; i<_thrusters.size(); i++)
         nonAeroWgt += ((ThrustRec*)_thrusters.get(i))->mass;
 
     // Rescale to the specified empty weight
     float wscale = (_emptyWeight-nonAeroWgt)/aeroWgt;
-    for(int i=firstMass; i<body->numMasses(); i++)
+    for(i=firstMass; i<body->numMasses(); i++)
         body->setMass(i, body->getMass(i)*wscale);
 
     // Add the thruster masses
-    for(int i=0; i<_thrusters.size(); i++) {
+    for(i=0; i<_thrusters.size(); i++) {
         ThrustRec* t = (ThrustRec*)_thrusters.get(i);
         body->addMass(t->mass, t->cg);
     }
 
     // Add the tanks, empty for now.
     float totalFuel = 0;
-    for(int i=0; i<_tanks.size(); i++) { 
+    for(i=0; i<_tanks.size(); i++) { 
         Tank* t = (Tank*)_tanks.get(i); 
         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();
 
     // Add surfaces for the landing gear.
-    for(int i=0; i<_gears.size(); i++)
+    for(i=0; i<_gears.size(); i++)
         compileGear((GearRec*)_gears.get(i));
 
     // The Thruster objects
-    for(int i=0; i<_thrusters.size(); i++) {
+    for(i=0; i<_thrusters.size(); i++) {
         ThrustRec* tr = (ThrustRec*)_thrusters.get(i);
         tr->handle = _model.addThruster(tr->thruster);
     }
@@ -492,13 +581,14 @@ void Airplane::compile()
     // Ground effect
     float gepos[3];
     float gespan = _wing->getGroundEffect(gepos);
-    _model.setGroundEffect(gepos, gespan, .3);
+    _model.setGroundEffect(gepos, gespan, 0.3f);
 
     solveGear();
     solve();
 
-    // 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()
@@ -512,30 +602,31 @@ void Airplane::solveGear()
     // "buffer" to keep things from blowing up with aircraft with a
     // single gear very near the c.g. (AV-8, for example).
     float total = 0;
-    for(int i=0; i<_gears.size(); i++) {
+    int i;
+    for(i=0; i<_gears.size(); i++) {
         GearRec* gr = (GearRec*)_gears.get(i);
         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;
     }
 
     // Renormalize so they sum to 1
-    for(int i=0; i<_gears.size(); i++)
+    for(i=0; i<_gears.size(); i++)
         ((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
     // a 3 degree approach.
-    float descentRate = 3*_approachSpeed/19.1;
+    float descentRate = 3.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(int i=0; i<_gears.size(); i++) {
+    for(i=0; i<_gears.size(); i++) {
         GearRec* gr = (GearRec*)_gears.get(i);
         float e = energy * gr->wgt;
         float comp[3];
@@ -551,14 +642,23 @@ void Airplane::solveGear()
         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->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();
     }
 }
 
 void Airplane::stabilizeThrust()
 {
-    for(int i=0; i<_thrusters.size(); i++)
+    int i;
+    for(i=0; i<_thrusters.size(); i++)
        _model.getThruster(i)->stabilize();
 }
 
@@ -570,32 +670,32 @@ void Airplane::runCruise()
 
     // The control configuration
     _controls.reset();
-    for(int i=0; i<_cruiseControls.size(); i++) {
+    int i;
+    for(i=0; i<_cruiseControls.size(); i++) {
        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);
    
     // Set up the thruster parameters and iterate until the thrust
     // stabilizes.
-    for(int i=0; i<_thrusters.size(); i++) {
+    for(i=0; i<_thrusters.size(); i++) {
        Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
        t->setWind(wind);
        t->setAir(_cruiseP, _cruiseT);
     }
     stabilizeThrust();
 
+    updateGearState();
+
     // Precompute thrust in the model, and calculate aerodynamic forces
     _model.getBody()->reset();
     _model.initIteration();
@@ -610,11 +710,12 @@ void Airplane::runApproach()
 
     // The control configuration
     _controls.reset();
-    for(int i=0; i<_approachControls.size(); i++) {
+    int i;
+    for(i=0; i<_approachControls.size(); i++) {
        Control* c = (Control*)_approachControls.get(i);
        _controls.setInput(c->control, c->val);
     }
-    _controls.applyControls();
+    _controls.applyControls(1000000);
 
     // The local wind
     float wind[3];
@@ -622,20 +723,19 @@ void Airplane::runApproach()
     Math::vmul33(_approachState.orient, wind, wind);
     
     // Approach is by convention at 20% tank capacity
-    setFuelFraction(0.2);
-
-    // Gear are down
-    setGearState(true, 100000);
+    setFuelFraction(0.2f);
 
     // Run the thrusters until they get to a stable setting.  FIXME:
     // this is lots of wasted work.
-    for(int i=0; i<_thrusters.size(); i++) {
+    for(i=0; i<_thrusters.size(); i++) {
        Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
        t->setWind(wind);
        t->setAir(_approachP, _approachT);
     }
     stabilizeThrust();
 
+    updateGearState();
+
     // Precompute thrust in the model, and calculate aerodynamic forces
     _model.getBody()->reset();
     _model.initIteration();
@@ -648,11 +748,12 @@ void Airplane::applyDragFactor(float factor)
     _dragFactor *= applied;
     _wing->setDragScale(_wing->getDragScale() * applied);
     _tail->setDragScale(_tail->getDragScale() * applied);
-    for(int i=0; i<_vstabs.size(); i++) {
+    int i;
+    for(i=0; i<_vstabs.size(); i++) {
        Wing* w = (Wing*)_vstabs.get(i);
        w->setDragScale(w->getDragScale() * applied);
     }
-    for(int i=0; i<_surfs.size(); i++) {
+    for(i=0; i<_surfs.size(); i++) {
        Surface* s = (Surface*)_surfs.get(i);
        s->setTotalDrag(s->getTotalDrag() * applied);
     }
@@ -664,7 +765,8 @@ void Airplane::applyLiftRatio(float factor)
     _liftRatio *= applied;
     _wing->setLiftRatio(_wing->getLiftRatio() * applied);
     _tail->setLiftRatio(_tail->getLiftRatio() * applied);
-    for(int i=0; i<_vstabs.size(); i++) {
+    int i;
+    for(i=0; i<_vstabs.size(); i++) {
         Wing* w = (Wing*)_vstabs.get(i);
         w->setLiftRatio(w->getLiftRatio() * applied);
     }
@@ -686,7 +788,7 @@ 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;
@@ -713,6 +815,9 @@ void Airplane::solve()
        // Run an approach iteration, and do likewise
        runApproach();
 
+       _model.getBody()->getAngularAccel(tmp);
+       float apitch0 = tmp[1];
+
        _model.getBody()->getAccel(tmp);
        float alift = _approachWeight * tmp[2];
 
@@ -733,7 +838,7 @@ void Airplane::solve()
        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);
@@ -749,7 +854,23 @@ void Airplane::solve()
             return;
         }
 
-       // And apply:
+        // 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.0001f;
+        _approachElevator.val += ELEVDIDDLE;
+        runApproach();
+        _approachElevator.val -= ELEVDIDDLE;
+
+       _model.getBody()->getAngularAccel(tmp);
+       float 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.
+
        applyDragFactor(dragFactor);
        applyLiftRatio(liftFactor);
 
@@ -760,15 +881,20 @@ void Airplane::solve()
            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 += 0.5f*aoaDelta;
+       _tailIncidence += 0.5f*tailDelta;
+        _approachElevator.val += 0.5f*elevDelta;
        
-       _cruiseAoA = clamp(_cruiseAoA, -.174, .174);
-       _tailIncidence = clamp(_tailIncidence, -.174, .174);
-
-        if(dragFactor < 1.00001 && liftFactor < 1.00001 &&
-           aoaDelta < .000017   && tailDelta < .000017)
+       _cruiseAoA = clamp(_cruiseAoA, -0.174f, 0.174f);
+       _tailIncidence = clamp(_tailIncidence, -0.174f, 0.174f);
+        _approachElevator.val = clamp(_approachElevator.val, -1.f, 1.f);
+
+        if(norm(dragFactor) < 1.00001 &&
+           norm(liftFactor) < 1.00001 &&
+           abs(aoaDelta) < .000017 &&
+           abs(tailDelta) < .000017 &&
+           abs(elevDelta) < 0.00001)
         {
             break;
         }