1 #include <Main/fg_props.hxx>
4 #include "Atmosphere.hpp"
7 // Copied from McCormick, who got it from "The ARDC Model Atmosphere"
8 // Note that there's an error in the text in the first entry,
9 // McCormick lists 299.16/101325/1.22500, but those don't agree with
10 // R=287. I chose to correct the temperature to 288.20, since 79F is
11 // pretty hot for a "standard" atmosphere.
12 // meters kelvin Pa kg/m^3
13 float Atmosphere::data[][4] = {{ 0.0f, 288.20f, 101325.0f, 1.22500f },
14 { 900.0f, 282.31f, 90971.0f, 1.12260f },
15 { 1800.0f, 276.46f, 81494.0f, 1.02690f },
16 { 2700.0f, 270.62f, 72835.0f, 0.93765f },
17 { 3600.0f, 264.77f, 64939.0f, 0.85445f },
18 { 4500.0f, 258.93f, 57752.0f, 0.77704f },
19 { 5400.0f, 253.09f, 51226.0f, 0.70513f },
20 { 6300.0f, 247.25f, 45311.0f, 0.63845f },
21 { 7200.0f, 241.41f, 39963.0f, 0.57671f },
22 { 8100.0f, 235.58f, 35140.0f, 0.51967f },
23 { 9000.0f, 229.74f, 30800.0f, 0.46706f },
24 { 9900.0f, 223.91f, 26906.0f, 0.41864f },
25 { 10800.0f, 218.08f, 23422.0f, 0.37417f },
26 { 11700.0f, 216.66f, 20335.0f, 0.32699f },
27 { 12600.0f, 216.66f, 17654.0f, 0.28388f },
28 { 13500.0f, 216.66f, 15327.0f, 0.24646f },
29 { 14400.0f, 216.66f, 13308.0f, 0.21399f },
30 { 15300.0f, 216.66f, 11555.0f, 0.18580f },
31 { 16200.0f, 216.66f, 10033.0f, 0.16133f },
32 { 17100.0f, 216.66f, 8712.0f, 0.14009f },
33 { 18000.0f, 216.66f, 7565.0f, 0.12165f },
34 { 18900.0f, 216.66f, 6570.0f, 0.10564f }};
36 // Universal gas constant for air, in SI units. P = R * rho * T.
37 // P in pascals (N/m^2), rho is kg/m^3, T in kelvin.
38 const float R = 287.1f;
40 // Specific heat ratio for air, at "low" temperatures.
41 const float GAMMA = 1.4f;
43 float Atmosphere::getStdTemperature(float alt)
45 if (fgGetBool("/environment/params/control-fdm-atmosphere"))
46 return fgGetDouble("/environment/temperature-degC") + 273.15;
48 return getRecord(alt, 1);
51 float Atmosphere::getStdPressure(float alt)
53 if (fgGetBool("/environment/params/control-fdm-atmosphere"))
54 return fgGetDouble("/environment/pressure-inhg") * 3386.39;
56 return getRecord(alt, 2);
59 float Atmosphere::getStdDensity(float alt)
61 if (fgGetBool("/environment/params/control-fdm-atmosphere"))
62 return fgGetDouble("/environment/density-slugft3") * 515.378;
64 return getRecord(alt, 3);
67 float Atmosphere::calcVEAS(float spd, float pressure, float temp)
69 static float rho0 = getStdDensity(0);
70 float densityRatio = calcDensity(pressure, temp) / rho0;
71 return spd * Math::sqrt(densityRatio);
74 float Atmosphere::calcVCAS(float spd, float pressure, float temp)
76 // Stolen shamelessly from JSBSim. Constants that appear:
78 // 5/12 == 1/(gamma+1)
82 // 144/25 == (gamma+1)^2
84 float m2 = calcMach(spd, temp);
87 float cp; // pressure coefficient
89 // (1+(mach^2)/5)^(gamma/(gamma-1))
90 cp = Math::pow(1+0.2*m2, 3.5);
92 float tmp0 = ((144.0f/25.0f) * m2) / (28.0f/5.0f*m2 - 4.0f/5.0f);
93 float tmp1 = ((14.0f/5.0f) * m2 - (2.0f/5.0f)) * (5.0f/12.0f);
94 cp = Math::pow(tmp0, 3.5) * tmp1;
97 // Conditions at sea level
98 float p0 = getStdPressure(0);
99 float rho0 = getStdDensity(0);
101 float tmp = Math::pow((pressure/p0)*(cp-1) + 1, (2/7.));
102 return Math::sqrt((7*p0/rho0)*(tmp-1));
105 float Atmosphere::calcDensity(float pressure, float temp)
107 return pressure / (R * temp);
110 float Atmosphere::calcMach(float spd, float temp)
112 return spd / Math::sqrt(GAMMA * R * temp);
115 void Atmosphere::calcStaticAir(float p0, float t0, float d0, float v,
116 float* pOut, float* tOut, float* dOut)
118 const static float C0 = ((GAMMA-1)/(2*R*GAMMA));
119 const static float C1 = 1/(GAMMA-1);
121 *tOut = t0 + (v*v) * C0;
122 *dOut = d0 * Math::pow(*tOut / t0, C1);
123 *pOut = (*dOut) * R * (*tOut);
126 float Atmosphere::getRecord(float alt, int recNum)
128 int hi = (sizeof(data) / (4*sizeof(float))) - 1;
131 // safety valve, clamp to the edges of the table
132 if(alt < data[0][0]) hi=1;
133 else if(alt > data[hi][0]) lo = hi-1;
137 if(hi-lo == 1) break;
138 int mid = (hi+lo)>>1;
139 if(alt < data[mid][0]) hi = mid;
144 float frac = (alt - data[lo][0])/(data[hi][0] - data[lo][0]);
145 float a = data[lo][recNum];
146 float b = data[hi][recNum];
147 return a + frac * (b-a);
150 }; // namespace yasim