src/FDM/YASim/Atmosphere.cpp

   1 #include "Math.hpp"
   2 #include "Atmosphere.hpp"
   3 namespace yasim {
   4
   5 // Copied from McCormick, who got it from "The ARDC Model Atmosphere"
   6 // Note that there's an error in the text in the first entry,
   7 // McCormick lists 299.16/101325/1.22500, but those don't agree with
   8 // R=287.  I chose to correct the temperature to 288.20, since 79F is
   9 // pretty hot for a "standard" atmosphere.
  10 //                             meters   kelvin     kPa   kg/m^3
  11 float Atmosphere::data[][4] = {{ 0,     288.20, 101325, 1.22500 },
  12                                { 900,   282.31,  90971, 1.12260 },
  13                                { 1800,  276.46,  81494, 1.02690 },
  14                                { 2700,  270.62,  72835, 0.93765 },
  15                                { 3600,  264.77,  64939, 0.85445 },
  16                                { 4500,  258.93,  57752, 0.77704 },
  17                                { 5400,  253.09,  51226, 0.70513 },
  18                                { 6300,  247.25,  45311, 0.63845 },
  19                                { 7200,  241.41,  39963, 0.57671 },
  20                                { 8100,  235.58,  35140, 0.51967 },
  21                                { 9000,  229.74,  30800, 0.46706 },
  22                                { 9900,  223.91,  26906, 0.41864 },
  23                                { 10800, 218.08,  23422, 0.37417 },
  24                                { 11700, 216.66,  20335, 0.32699 },
  25                                { 12600, 216.66,  17654, 0.28388 },
  26                                { 13500, 216.66,  15327, 0.24646 },
  27                                { 14400, 216.66,  13308, 0.21399 },
  28                                { 15300, 216.66,  11555, 0.18580 },
  29                                { 16200, 216.66,  10033, 0.16133 },
  30                                { 17100, 216.66,   8712, 0.14009 },
  31                                { 18000, 216.66,   7565, 0.12165 },
  32                                { 18900, 216.66,   6570, 0.10564 }};
  33
  34 float Atmosphere::getStdTemperature(float alt)
  35 {
  36     return getRecord(alt, 1);
  37 }
  38
  39 float Atmosphere::getStdPressure(float alt)
  40 {
  41     return getRecord(alt, 2);
  42 }
  43
  44 float Atmosphere::getStdDensity(float alt)
  45 {
  46     return getRecord(alt, 3);
  47 }
  48
  49 float Atmosphere::calcVEAS(float spd, float pressure, float temp)
  50 {
  51     return 0; //FIXME
  52 }
  53
  54 float Atmosphere::calcVCAS(float spd, float pressure, float temp)
  55 {
  56     // Stolen shamelessly from JSBSim.  Constants that appear:
  57     //   2/5  == gamma-1
  58     //   5/12 == 1/(gamma+1)
  59     //   4/5  == 2*(gamma-1)
  60     //  14/5  == 2*gamma
  61     //  28/5  == 4*gamma
  62     // 144/25 == (gamma+1)^2
  63
  64     float m2 = calcMach(spd, temp);
  65     m2 = m2*m2; // mach^2
  66
  67     float cp; // pressure coefficient
  68     if(m2 < 1) {
  69         // (1+(mach^2)/5)^(gamma/(gamma-1))
  70         cp = Math::pow(1+0.2*m2, 3.5);
  71     } else {
  72         float tmp0 = ((144/25.) * m2) / (28/5.*m2 - 4/5.);
  73         float tmp1 = ((14/5.) * m2 - (2/5.)) * (5/12.);
  74         cp = Math::pow(tmp0, 3.5) * tmp1;
  75     }
  76
  77     // Conditions at sea level
  78     float p0 = getStdPressure(0);
  79     float rho0 = getStdDensity(0);
  80
  81     float tmp = Math::pow((pressure/p0)*(cp-1) + 1, (2/7.));
  82     return Math::sqrt((7*p0/rho0)*(tmp-1));
  83 }
  84
  85 float Atmosphere::calcDensity(float pressure, float temp)
  86 {
  87     // P = rho*R*T, R == 287 kPa*m^3 per kg*kelvin for air
  88     return pressure / (287 * temp);
  89 }
  90
  91 float Atmosphere::calcMach(float spd, float temp)
  92 {
  93     return spd / Math::sqrt(1.4 * 287 * temp);
  94 }
  95
  96 float Atmosphere::getRecord(float alt, int recNum)
  97 {
  98     int hi = (sizeof(data) / (4*sizeof(float))) - 1;
  99     int lo = 0;
 100
 101     // safety valve, clamp to the edges of the table
 102     if(alt < data[0][0])       hi=1;
 103     else if(alt > data[hi][0]) lo = hi-1;
 104
 105     // binary search
 106     while(1) {
 107         if(hi-lo == 1) break;
 108         int mid = (hi+lo)>>1;
 109         if(alt < data[mid][0]) hi = mid;
 110         else lo = mid;
 111     }
 112
 113     // interpolate
 114     float frac = (alt - data[lo][0])/(data[hi][0] - data[lo][0]);
 115     float a = data[lo][recNum];
 116     float b = data[hi][recNum];
 117     return a + frac * (b-a);
 118 }
 119
 120 }; // namespace yasim