src/FDM/YASim/yasim-test.cpp

   1 #include <stdio.h>
   2
   3 #include <cstring>
   4 #include <cstdlib>
   5
   6 #include <simgear/props/props.hxx>
   7 #include <simgear/xml/easyxml.hxx>
   8
   9 #include "FGFDM.hpp"
  10 #include "Atmosphere.hpp"
  11 #include "Airplane.hpp"
  12
  13 #include <simgear/math/SGMath.hxx>
  14
  15 using namespace yasim;
  16
  17 // Stubs.  Not needed by a batch program, but required to link.
  18 bool fgSetFloat (const char * name, float val) { return false; }
  19 bool fgSetBool(char const * name, bool val) { return false; }
  20 bool fgGetBool(char const * name, bool def) { return false; }
  21 bool fgSetString(char const * name, char const * str) { return false; }
  22 SGPropertyNode* fgGetNode (const char * path, bool create) { return 0; }
  23 SGPropertyNode* fgGetNode (const char * path, int i, bool create) { return 0; }
  24 float fgGetFloat (const char * name, float defaultValue) { return 0; }
  25 double fgGetDouble (const char * name, double defaultValue = 0.0) { return 0; }
  26 bool fgSetDouble (const char * name, double defaultValue = 0.0) { return 0; }
  27
  28 static const float RAD2DEG = 57.2957795131;
  29 static const float DEG2RAD = 0.0174532925199;
  30 static const float KTS2MPS = 0.514444444444;
  31
  32
  33 // Generate a graph of lift, drag and L/D against AoA at the specified
  34 // speed and altitude.  The result is a space-separated file of
  35 // numbers: "aoa lift drag LD" (aoa in degrees, lift and drag in
  36 // G's).  You can use this in gnuplot like so (assuming the output is
  37 // in a file named "dat":
  38 /*
  39  plot "dat" using 1:2 with lines title 'lift', \
  40       "dat" using 1:3 with lines title 'drag', \
  41       "dat" using 1:4 with lines title 'LD'
  42 */
  43 void yasim_graph(Airplane* a, float alt, float kts)
  44 {
  45     Model* m = a->getModel();
  46     State s;
  47
  48     m->setAir(Atmosphere::getStdPressure(alt),
  49               Atmosphere::getStdTemperature(alt),
  50               Atmosphere::getStdDensity(alt));
  51     m->getBody()->recalc();
  52
  53     for(int deg=-179; deg<=179; deg++) {
  54         float aoa = deg * DEG2RAD;
  55         Airplane::setupState(aoa, kts * KTS2MPS, 0 ,&s);
  56         m->getBody()->reset();
  57         m->initIteration();
  58         m->calcForces(&s);
  59
  60         float acc[3];
  61         m->getBody()->getAccel(acc);
  62         Math::tmul33(s.orient, acc, acc);
  63
  64         float drag = acc[0] * (-1/9.8);
  65         float lift = 1 + acc[2] * (1/9.8);
  66
  67         printf("%d %g %g %g\n", deg, lift, drag, lift/drag);
  68     }
  69 }
  70
  71 int usage()
  72 {
  73     fprintf(stderr, "Usage: yasim <ac.xml> [-g [-a alt] [-s kts]]\n");
  74     return 1;
  75 }
  76
  77 int main(int argc, char** argv)
  78 {
  79     FGFDM* fdm = new FGFDM();
  80     Airplane* a = fdm->getAirplane();
  81
  82     if(argc < 2) return usage();
  83
  84     // Read
  85     try {
  86         string file = argv[1];
  87         readXML(file, *fdm);
  88     } catch (const sg_exception &e) {
  89         printf("XML parse error: %s (%s)\n",
  90                e.getFormattedMessage().c_str(), e.getOrigin());
  91     }
  92
  93     // ... and run
  94     a->compile();
  95     if(a->getFailureMsg())
  96         printf("SOLUTION FAILURE: %s\n", a->getFailureMsg());
  97
  98     if(!a->getFailureMsg() && argc > 2 && strcmp(argv[2], "-g") == 0) {
  99         float alt = 5000, kts = 100;
 100         for(int i=3; i<argc; i++) {
 101             if     (std::strcmp(argv[i], "-a") == 0) alt = std::atof(argv[++i]);
 102             else if(std::strcmp(argv[i], "-s") == 0) kts = std::atof(argv[++i]);
 103             else return usage();
 104         }
 105         yasim_graph(a, alt, kts);
 106     } else {
 107         float aoa = a->getCruiseAoA() * RAD2DEG;
 108         float tail = -1 * a->getTailIncidence() * RAD2DEG;
 109         float drag = 1000 * a->getDragCoefficient();
 110         float cg[3];
 111         a->getModel()->getBody()->getCG(cg);
 112         a->getModel()->getBody()->recalc();
 113
 114         float SI_inertia[9];
 115         a->getModel()->getBody()->getInertiaMatrix(SI_inertia);
 116
 117         printf("Solution results:");
 118         printf("       Iterations: %d\n", a->getSolutionIterations());
 119         printf(" Drag Coefficient: %f\n", drag);
 120         printf("       Lift Ratio: %f\n", a->getLiftRatio());
 121         printf("       Cruise AoA: %f\n", aoa);
 122         printf("   Tail Incidence: %f\n", tail);
 123         printf("Approach Elevator: %f\n", a->getApproachElevator());
 124         printf("               CG: x:%.3f, y:%.3f, z:%.3f\n\n", cg[0], cg[1], cg[2]);
 125         printf("  Inertia tensor : %.3f, %.3f, %.3f\n", SI_inertia[0], SI_inertia[1], SI_inertia[2]);
 126         printf("        [kg*m^2]   %.3f, %.3f, %.3f\n", SI_inertia[3], SI_inertia[4], SI_inertia[5]);
 127         printf("     Origo at CG   %.3f, %.3f, %.3f\n", SI_inertia[6], SI_inertia[7], SI_inertia[8]);
 128     }
 129     delete fdm;
 130     return 0;
 131 }