Add a lift/drag vs. AoA graph option to the yasim tool, which

generates a data file of aerodynamic lift and drag (and L/D) against
AoA at a specified speed and altitude through a full circle.  Wrote it
to track down the YF-23 superthrust issue, but it wasn't any help.
All the forces look fine.

src/FDM/YASim/Airplane.hpp

@@ -90,6 +90,8 @@ public:
     float getApproachElevator() { return _approachElevator.val; }
     char* getFailureMsg();
 
+    static void setupState(float aoa, float speed, State* s); // utility
+
 private:
     struct Tank { float pos[3]; float cap; float fill;
 	          float density; int handle; };
@@ -104,7 +106,6 @@ private:
 
     void runCruise();
     void runApproach();
-    void setupState(float aoa, float speed, State* s);
     void solveGear();
     void solve();
     void solveHelicopter();

src/FDM/YASim/yasim-test.cpp

@@ -4,6 +4,7 @@
 #include <simgear/xml/easyxml.hxx>
 
 #include "FGFDM.hpp"
+#include "Atmosphere.hpp"
 #include "Airplane.hpp"
 
 using namespace yasim;
@@ -19,12 +20,61 @@ float fgGetDouble (const char * name, double defaultValue) { return 0; }
 float fgSetDouble (const char * name, double defaultValue) { return 0; }
 
 static const float RAD2DEG = 57.2957795131;
+static const float DEG2RAD = 0.0174532925199;
+static const float KTS2MPS = 0.514444444444;
+
+
+// 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
+// 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":
+//
+// plot "dat" using 1:2 with lines title 'lift', \
+//      "dat" using 1:3 with lines title 'drag', \
+//      "dat" using 1:4 with lines title 'LD'
+//
+void yasim_graph(Airplane* a, float alt, float kts)
+{
+    Model* m = a->getModel();
+    State s;
+
+    m->setAir(Atmosphere::getStdPressure(alt),
+              Atmosphere::getStdTemperature(alt),
+              Atmosphere::getStdDensity(alt));
+    m->getBody()->recalc();
+
+    for(int deg=-179; deg<=179; deg++) {
+        float aoa = deg * DEG2RAD;
+        Airplane::setupState(aoa, kts * KTS2MPS, &s);
+        m->getBody()->reset();
+        m->initIteration();
+        m->calcForces(&s);
+
+        float acc[3];
+	m->getBody()->getAccel(acc);
+        Math::tmul33(s.orient, acc, acc);
+
+        float drag = acc[0] * (-1/9.8);
+        float lift = 1 + acc[2] * (1/9.8);
+
+        printf("%d %g %g %g\n", deg, lift, drag, lift/drag);
+    }
+}
+
+int usage()
+{
+    fprintf(stderr, "Usage: yasim <ac.xml> [-g [-a alt] [-s kts]]\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];
@@ -36,26 +86,33 @@ int main(int argc, char** argv)
 
     // ... and run
     a->compile();
-
-    float aoa = a->getCruiseAoA() * RAD2DEG;
-    float tail = -1 * a->getTailIncidence() * RAD2DEG;
-    float drag = 1000 * a->getDragCoefficient();
-    float cg[3];
-    a->getModel()->getBody()->getCG(cg);
-
-    printf("Solution results:");
-    printf("       Iterations: %d\n", a->getSolutionIterations());
-    printf(" Drag Coefficient: %f\n", drag);
-    printf("       Lift Ratio: %f\n", a->getLiftRatio());
-    printf("       Cruise AoA: %f\n", aoa);
-    printf("   Tail Incidence: %f\n", tail);
-    printf("Approach Elevator: %f\n", a->getApproachElevator());
-    printf("               CG: %.3f, %.3f, %.3f\n", cg[0], cg[1], cg[2]);
-
     if(a->getFailureMsg())
         printf("SOLUTION FAILURE: %s\n", a->getFailureMsg());
 
+    if(!a->getFailureMsg() && argc > 2 && strcmp(argv[2], "-g") == 0) {
+        float alt = 5000, kts = 100;
+        for(int i=3; i<argc; i++) {
+            if     (strcmp(argv[i], "-a") == 0) alt = atof(argv[++i]);
+            else if(strcmp(argv[i], "-s") == 0) kts = atof(argv[++i]);
+            else return usage();
+        }
+        yasim_graph(a, alt, kts);
+    } else {
+        float aoa = a->getCruiseAoA() * RAD2DEG;
+        float tail = -1 * a->getTailIncidence() * RAD2DEG;
+        float drag = 1000 * a->getDragCoefficient();
+        float cg[3];
+        a->getModel()->getBody()->getCG(cg);
+        
+        printf("Solution results:");
+        printf("       Iterations: %d\n", a->getSolutionIterations());
+        printf(" Drag Coefficient: %f\n", drag);
+        printf("       Lift Ratio: %f\n", a->getLiftRatio());
+        printf("       Cruise AoA: %f\n", aoa);
+        printf("   Tail Incidence: %f\n", tail);
+        printf("Approach Elevator: %f\n", a->getApproachElevator());
+        printf("               CG: %.3f, %.3f, %.3f\n", cg[0], cg[1], cg[2]);
+    }
     delete fdm;
-
     return 0;
 }