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      Pa   kg/m^3
  11 float Atmosphere::data[][4] = {{ 0.0f,     288.20f, 101325.0f, 1.22500f },
  12                                {   900.0f, 282.31f,  90971.0f, 1.12260f },
  13                                {  1800.0f, 276.46f,  81494.0f, 1.02690f },
  14                                {  2700.0f, 270.62f,  72835.0f, 0.93765f },
  15                                {  3600.0f, 264.77f,  64939.0f, 0.85445f },
  16                                {  4500.0f, 258.93f,  57752.0f, 0.77704f },
  17                                {  5400.0f, 253.09f,  51226.0f, 0.70513f },
  18                                {  6300.0f, 247.25f,  45311.0f, 0.63845f },
  19                                {  7200.0f, 241.41f,  39963.0f, 0.57671f },
  20                                {  8100.0f, 235.58f,  35140.0f, 0.51967f },
  21                                {  9000.0f, 229.74f,  30800.0f, 0.46706f },
  22                                {  9900.0f, 223.91f,  26906.0f, 0.41864f },
  23                                { 10800.0f, 218.08f,  23422.0f, 0.37417f },
  24                                { 11700.0f, 216.66f,  20335.0f, 0.32699f },
  25                                { 12600.0f, 216.66f,  17654.0f, 0.28388f },
  26                                { 13500.0f, 216.66f,  15327.0f, 0.24646f },
  27                                { 14400.0f, 216.66f,  13308.0f, 0.21399f },
  28                                { 15300.0f, 216.66f,  11555.0f, 0.18580f },
  29                                { 16200.0f, 216.66f,  10033.0f, 0.16133f },
  30                                { 17100.0f, 216.66f,   8712.0f, 0.14009f },
  31                                { 18000.0f, 216.66f,   7565.0f, 0.12165f },
  32                                { 18900.0f, 216.66f,   6570.0f, 0.10564f }};
  33
  34 // Universal gas constant for air, in SI units.  P = R * rho * T.
  35 // P in pascals (N/m^2), rho is kg/m^3, T in kelvin.
  36 const float R = 287.1f;
  37
  38 // Specific heat ratio for air, at "low" temperatures.
  39 const float GAMMA = 1.4f;
  40
  41 float Atmosphere::getStdTemperature(float alt)
  42 {
  43     return getRecord(alt, 1);
  44 }
  45
  46 float Atmosphere::getStdPressure(float alt)
  47 {
  48     return getRecord(alt, 2);
  49 }
  50
  51 float Atmosphere::getStdDensity(float alt)
  52 {
  53     return getRecord(alt, 3);
  54 }
  55
  56 float Atmosphere::calcVEAS(float spd,
  57                            float pressure, float temp, float density)
  58 {
  59     static float rho0 = getStdDensity(0);
  60     float densityRatio = density / rho0;
  61     return spd * Math::sqrt(densityRatio);
  62 }
  63
  64 float Atmosphere::calcVCAS(float spd, float pressure, float temp)
  65 {
  66     // Stolen shamelessly from JSBSim.  Constants that appear:
  67     //   2/5  == gamma-1
  68     //   5/12 == 1/(gamma+1)
  69     //   4/5  == 2*(gamma-1)
  70     //  14/5  == 2*gamma
  71     //  28/5  == 4*gamma
  72     // 144/25 == (gamma+1)^2
  73
  74     float m2 = calcMach(spd, temp);
  75     m2 = m2*m2; // mach^2
  76
  77     float cp; // pressure coefficient
  78     if(m2 < 1) {
  79         // (1+(mach^2)/5)^(gamma/(gamma-1))
  80         cp = Math::pow(1+0.2*m2, 3.5);
  81     } else {
  82         float tmp0 = ((144.0f/25.0f) * m2) / (28.0f/5.0f*m2 - 4.0f/5.0f);
  83         float tmp1 = ((14.0f/5.0f) * m2 - (2.0f/5.0f)) * (5.0f/12.0f);
  84         cp = Math::pow(tmp0, 3.5) * tmp1;
  85     }
  86
  87     // Conditions at sea level
  88     float p0 = getStdPressure(0);
  89     float rho0 = getStdDensity(0);
  90
  91     float tmp = Math::pow((pressure/p0)*(cp-1) + 1, (2/7.));
  92     return Math::sqrt((7*p0/rho0)*(tmp-1));
  93 }
  94
  95 float Atmosphere::calcStdDensity(float pressure, float temp)
  96 {
  97     return pressure / (R * temp);
  98 }
  99
 100 float Atmosphere::calcMach(float spd, float temp)
 101 {
 102     return spd / Math::sqrt(GAMMA * R * temp);
 103 }
 104
 105 float Atmosphere::spdFromMach(float mach, float temp)
 106 {
 107     return mach * Math::sqrt(GAMMA * R * temp);
 108 }
 109
 110 float Atmosphere::spdFromVCAS(float vcas, float pressure, float temp)
 111 {
 112                                 // FIXME: does not account for supersonic
 113     float p0 = getStdPressure(0);
 114     float rho0 = getStdDensity(0);
 115
 116     float tmp = (vcas*vcas)/(7*p0/rho0) + 1;
 117     float cp = ((Math::pow(tmp,(7/2.))-1)/(pressure/p0)) + 1;
 118
 119     float m2 = (Math::pow(cp,(1/3.5))-1)/0.2;
 120     float vtas= spdFromMach(Math::sqrt(m2), temp);
 121     return vtas;
 122 }
 123
 124 void Atmosphere::calcStaticAir(float p0, float t0, float d0, float v,
 125                                float* pOut, float* tOut, float* dOut)
 126 {
 127     const static float C0 = ((GAMMA-1)/(2*R*GAMMA));
 128     const static float C1 = 1/(GAMMA-1);
 129
 130     *tOut = t0 + (v*v) * C0;
 131     *dOut = d0 * Math::pow(*tOut / t0, C1);
 132     *pOut = (*dOut) * R * (*tOut);
 133 }
 134
 135 float Atmosphere::getRecord(float alt, int recNum)
 136 {
 137     int hi = (sizeof(data) / (4*sizeof(float))) - 1;
 138     int lo = 0;
 139
 140     // safety valve, clamp to the edges of the table
 141     if(alt < data[0][0])       hi=1;
 142     else if(alt > data[hi][0]) lo = hi-1;
 143
 144     // binary search
 145     while(1) {
 146         if(hi-lo == 1) break;
 147         int mid = (hi+lo)>>1;
 148         if(alt < data[mid][0]) hi = mid;
 149         else lo = mid;
 150     }
 151
 152     // interpolate
 153     float frac = (alt - data[lo][0])/(data[hi][0] - data[lo][0]);
 154     float a = data[lo][recNum];
 155     float b = data[hi][recNum];
 156     return a + frac * (b-a);
 157 }
 158
 159 }; // namespace yasim