2 #include "Atmosphere.hpp"
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, 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 }};
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.1;
38 // Specific heat ratio for air, at "low" temperatures.
39 const float GAMMA = 1.4;
41 float Atmosphere::getStdTemperature(float alt)
43 return getRecord(alt, 1);
46 float Atmosphere::getStdPressure(float alt)
48 return getRecord(alt, 2);
51 float Atmosphere::getStdDensity(float alt)
53 return getRecord(alt, 3);
56 float Atmosphere::calcVEAS(float spd, float pressure, float temp)
58 static float rho0 = getStdDensity(0);
59 float densityRatio = calcDensity(pressure, temp) / rho0;
60 return spd * Math::sqrt(densityRatio);
63 float Atmosphere::calcVCAS(float spd, float pressure, float temp)
65 // Stolen shamelessly from JSBSim. Constants that appear:
67 // 5/12 == 1/(gamma+1)
71 // 144/25 == (gamma+1)^2
73 float m2 = calcMach(spd, temp);
76 float cp; // pressure coefficient
78 // (1+(mach^2)/5)^(gamma/(gamma-1))
79 cp = Math::pow(1+0.2*m2, 3.5);
81 float tmp0 = ((144/25.) * m2) / (28/5.*m2 - 4/5.);
82 float tmp1 = ((14/5.) * m2 - (2/5.)) * (5/12.);
83 cp = Math::pow(tmp0, 3.5) * tmp1;
86 // Conditions at sea level
87 float p0 = getStdPressure(0);
88 float rho0 = getStdDensity(0);
90 float tmp = Math::pow((pressure/p0)*(cp-1) + 1, (2/7.));
91 return Math::sqrt((7*p0/rho0)*(tmp-1));
94 float Atmosphere::calcDensity(float pressure, float temp)
96 return pressure / (R * temp);
99 float Atmosphere::calcMach(float spd, float temp)
101 return spd / Math::sqrt(GAMMA * R * temp);
104 void Atmosphere::calcStaticAir(float p0, float t0, float d0, float v,
105 float* pOut, float* tOut, float* dOut)
107 const static float C0 = ((GAMMA-1)/(2*R*GAMMA));
108 const static float C1 = 1/(GAMMA-1);
110 *tOut = t0 + (v*v) * C0;
111 *dOut = d0 * Math::pow(*tOut / t0, C1);
112 *pOut = (*dOut) * R * (*tOut);
115 float Atmosphere::getRecord(float alt, int recNum)
117 int hi = (sizeof(data) / (4*sizeof(float))) - 1;
120 // safety valve, clamp to the edges of the table
121 if(alt < data[0][0]) hi=1;
122 else if(alt > data[hi][0]) lo = hi-1;
126 if(hi-lo == 1) break;
127 int mid = (hi+lo)>>1;
128 if(alt < data[mid][0]) hi = mid;
133 float frac = (alt - data[lo][0])/(data[hi][0] - data[lo][0]);
134 float a = data[lo][recNum];
135 float b = data[hi][recNum];
136 return a + frac * (b-a);
139 }; // namespace yasim