1 #include <simgear/math/SGMath.hxx>
2 #include <simgear/debug/logstream.hxx>
4 #include "atmosphere.hxx"
9 const ISA_layer ISA_def[] = {
11 // id (m) (ft) (Pa) (inHg) (K) (C) (K/m) (K/ft)
12 ISA_layer(0, 0, 0, 101325, 29.92126, 288.15, 15.00, 0.0065, 0.0019812),
13 ISA_layer(1, 11000, 36089, 22632.1, 6.683246, 216.65, -56.50, 0, 0),
14 ISA_layer(2, 20000, 65616, 5474.89, 1.616734, 216.65, -56.50, -0.0010, -0.0003048),
15 ISA_layer(3, 32000, 104986, 868.019, 0.256326, 228.65, -44.50, -0.0028, -0.0008534),
16 ISA_layer(4, 47000, 154199, 110.906, 0.0327506, 270.65, -2.50, 0, 0),
17 ISA_layer(5, 51000, 167322, 66.9389, 0.0197670, 270.65, -2.50, 0.0028, 0.0008534),
18 ISA_layer(6, 71000, 232939, 3.95642, 0.00116833, 214.65, -58.50, 0.0020, 0.0006096),
19 ISA_layer(7, 80000, 262467, 0.88628, 0.000261718, 196.65, -76.50),
22 const int ISA_def_size(sizeof(ISA_def) / sizeof(ISA_layer));
25 // Pressure within a layer, as a function of height
26 // and of layer parameters.
27 // Physics model: standard or nonstandard atmosphere,
28 // depending on what parameters you pass in.
29 // Heights in meters, pressures in Pa.
30 // As always, $lambda is positive in the troposphere,
31 // and zero in the first part of the stratosphere.
38 const double lambda) {
39 using namespace atmodel;
41 // (1/N) &:=& {g\,m_M \over \lambda\,R}
42 double N = lambda * Rgas / mm / g;
43 return Pb * pow((Tb - lambda*(hh - hb)) / Tb, 1/N);
45 return Pb * exp(-g * mm / Rgas / Tb * (hh - hb));
49 // Pressure as a function of height.
50 // Valid below 32000 meters,
51 // i.e. troposphere and first two layers of stratosphere.
52 // Does not depend on any caching; can be used to
53 // *construct* caches and interpolation tables.
55 // Height in meters, pressure in pascals.
57 double old_p_vs_a(const double height) {
58 using namespace atmodel;
59 if (height <= 11000.) {
60 return P_layer(height, 0.0, ISA::P0, ISA::T0, ISA::lam0);
61 } else if (height <= 20000.) {
62 return P_layer(height, 11000., 22632.06, 216.65, 0.0);
63 } else if (height <= 32000.) {
64 return P_layer(height, 20000., 5474.89, 216.65, -0.001);
72 // Pressure within a layer, as a function of height.
73 // Physics model: standard or nonstandard atmosphere,
74 // depending on what parameters you pass in.
75 // Height in meters, pressures in pascals.
76 // As always, lapse is positive in the troposphere,
77 // and zero in the first part of the stratosphere.
79 double P_layer(const double height, const double href,
80 const double Pref, const double Tref,
82 using namespace atmodel;
84 double N = lapse * Rgas / mm / g;
85 return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
87 return Pref * exp(-g * mm / Rgas / Tref * (height - href));
92 // Temperature within a layer, as a function of height.
93 // Physics model: standard or nonstandard atmosphere
94 // depending on what parameters you pass in.
95 // $hh in meters, pressures in Pa.
96 // As always, $lambda is positive in the troposphere,
97 // and zero in the first part of the stratosphere.
103 const double lambda) {
104 return Tb - lambda*(hh - hb);
107 // Pressure and temperature as a function of height, Psl, and Tsl.
108 // heights in meters, pressures in Pa.
109 // Daisy chain version.
110 // We need "seed" values for sea-level pressure and temperature.
111 // In addition, for every layer, we need three things
112 // from the table: the reference height in that layer,
113 // the lapse in that layer, and the cap (if any) for that layer
114 // (which we take from the /next/ row of the table, if any).
115 tuple<double,double> PT_vs_hpt(
122 double hgt = ISA_def[0].height;
126 cout << "PT_vs_hpt: " << hh << " " << p0 << " " << t0 << endl;
130 for (const ISA_layer* pp = ISA_def; pp->lapse != -1; pp++, ii++) {
132 cout << "PT_vs_hpt: " << ii
133 << " height: " << pp->height
134 << " temp: " << pp->temp
135 << " lapse: " << pp->lapse
139 double lapse = pp->lapse;
140 // Stratosphere starts at a definite temperature,
141 // not a definite height:
143 xhgt = hgt + (t0 - (pp+1)->temp) / lapse;
144 } else if ((pp+1)->lapse != -1) {
145 xhgt = (pp+1)->height;
148 return make_tuple(P_layer(hh, hgt, p0, t0, lapse),
149 T_layer(hh, hgt, p0, t0, lapse));
151 p0 = P_layer(xhgt, hgt, p0, t0, lapse);
152 t0 = t0 - lapse * (xhgt - hgt);
156 // Should never get here.
157 SG_LOG(SG_GENERAL, SG_ALERT, "PT_vs_hpt: ran out of layers");
158 return make_tuple(d0,d0);
162 FGAtmoCache::FGAtmoCache() :
166 FGAtmoCache::~FGAtmoCache() {
172 // The following two routines are called "fake" because they
173 // bypass the exceedingly complicated layer model implied by
174 // the "weather conditioins" popup menu.
175 // For now we must bypass it for several reasons, including
176 // the fact that we don't have an "environment" object for
177 // the airport (only for the airplane).
178 // degrees C, height in feet
179 double FGAtmo::fake_T_vs_a_us(const double h_ft,
180 const double Tsl) const {
181 using namespace atmodel;
182 return Tsl - ISA::lam0 * h_ft * foot;
185 // Dewpoint. degrees C or K, height in feet
186 double FGAtmo::fake_dp_vs_a_us(const double dpsl, const double h_ft) {
187 const double dp_lapse(0.002); // [K/m] approximate
188 // Reference: http://en.wikipedia.org/wiki/Lapse_rate
189 return dpsl - dp_lapse * h_ft * atmodel::foot;
192 // Height as a function of pressure.
193 // Valid in the troposphere only.
194 double FGAtmo::a_vs_p(const double press, const double qnh) {
195 using namespace atmodel;
197 double nn = lam0 * Rgas / g / mm;
198 return T0 * ( pow(qnh/P0,nn) - pow(press/P0,nn) ) / lam0;
201 // force retabulation
202 void FGAtmoCache::tabulate() {
203 using namespace atmodel;
205 a_tvs_p = new SGInterpTable;
207 for (double hgt = -1000; hgt <= 32000;) {
209 make_tuple(ref(press), ref(temp)) = PT_vs_hpt(hgt);
210 a_tvs_p->addEntry(press / inHg, hgt / foot);
212 #ifdef DEBUG_EXPORT_P_H
214 char* fmt = " { %9.2f , %5.0f },";
215 if (press < 10000) fmt = " { %9.3f , %5.0f },";
216 snprintf(buf, 100, fmt, press, hgt);
227 // make sure cache is valid
228 void FGAtmoCache::cache() {
233 // Check the basic function,
234 // then compare against the interpolator.
235 void FGAtmoCache::check_model() {
254 for (int i = 0; ; i++) {
255 double height = hgts[i];
258 using namespace atmodel;
261 make_tuple(ref(press), ref(temp)) = PT_vs_hpt(height);
262 cout << "Height: " << height
263 << " \tpressure: " << press << endl;
265 << a_tvs_p->interpolate(press / inHg)*foot << endl;
269 //////////////////////////////////////////////////////////////////////
271 FGAltimeter::FGAltimeter() : kset(atmodel::ISA::P0), kft(0)
276 double FGAltimeter::reading_ft(const double p_inHg, const double set_inHg) {
277 using namespace atmodel;
278 double press_alt = a_tvs_p->interpolate(p_inHg);
279 double kollsman_shift = a_tvs_p->interpolate(set_inHg);
280 return (press_alt - kollsman_shift);
283 // Altimeter setting _in pascals_
284 // ... caller gets to convert to inHg or millibars
285 // Field elevation in m
286 // Field pressure in pascals
287 // Valid for fields within the troposphere only.
288 double FGAtmo::QNH(const double field_elev, const double field_press) {
289 using namespace atmodel;
291 // Equation derived in altimetry.htm
292 // exponent in QNH equation:
293 double nn = ISA::lam0 * Rgas / g / mm;
294 // pressure ratio factor:
295 double prat = pow(ISA::P0 / field_press, nn);
296 double rslt = field_press
297 * pow(1. + ISA::lam0 * field_elev / ISA::T0 * prat, 1./nn);
299 SG_LOG(SG_GENERAL, SG_ALERT, "QNH: elev: " << field_elev
300 << " press: " << field_press
304 << " rslt/inHG: " << rslt/inHg);
309 void FGAltimeter::dump_stack1(const double Tref) {
310 using namespace atmodel;
313 double Psl = P_layer(0, 0, ISA::P0, Tref, ISA::lam0);
314 snprintf(buf, bs, "Tref: %6.2f Psl: %5.0f = %7.4f",
315 Tref, Psl, Psl / inHg);
319 " %6s %6s %6s %6s %6s %6s %6s",
320 "A", "Aind", "Apr", "Aprind", "P", "Psl", "Qnh");
323 double hgts[] = {0, 2500, 5000, 7500, 10000, -9e99};
324 for (int ii = 0; ; ii++) {
325 double hgt_ft = hgts[ii];
326 double hgt = hgt_ft * foot;
329 double press = P_layer(hgt, 0, ISA::P0, Tref, ISA::lam0);
330 double qnhx = QNH(hgt, press) / inHg;
331 double qnh2 = SGMiscd::round(qnhx*100)/100;
333 double p_inHg = press / inHg;
334 double Aprind = reading_ft(p_inHg);
335 double Apr = a_vs_p(p_inHg*inHg) / foot;
336 double hind = reading_ft(p_inHg, qnh2);
338 " %6.0f %6.0f %6.0f %6.0f %6.2f %6.2f %6.2f",
339 hgt_ft, hind, Apr, Aprind, p_inHg, Psl/inHg, qnh2);
345 void FGAltimeter::dump_stack() {
346 using namespace atmodel;
347 cout << "........." << endl;
348 cout << "Size: " << sizeof(FGAtmo) << endl;
349 dump_stack1(ISA::T0);
350 dump_stack1(ISA::T0 - 20);