1 #include <simgear/math/SGMath.hxx>
3 #include "atmosphere.hxx"
8 FGAtmoCache::FGAtmoCache() :
12 FGAtmoCache::~FGAtmoCache() {
16 // Pressure as a function of height.
17 // Valid below 32000 meters,
18 // i.e. troposphere and first two layers of stratosphere.
19 // Does not depend on any caching; can be used to
20 // *construct* caches and interpolation tables.
22 // Height in meters, pressure in pascals.
24 double FGAtmo::p_vs_a(const double height) {
25 using namespace atmodel;
26 if (height <= 11000.) {
27 return P_layer(height, 0.0, ISA::P0, ISA::T0, ISA::lam0);
28 } else if (height <= 20000.) {
29 return P_layer(height, 11000., 22632.06, 216.65, 0.0);
30 } else if (height <= 32000.) {
31 return P_layer(height, 20000., 5474.89, 216.65, -0.001);
36 // degrees C, height in feet
37 double FGAtmo::fake_t_vs_a_us(const double h_ft) {
38 using namespace atmodel;
39 return ISA::T0 - ISA::lam0 * h_ft * foot - freezing;
42 // Dewpoint. degrees C or K, height in feet
43 double FGAtmo::fake_dp_vs_a_us(const double dpsl, const double h_ft) {
44 const double dp_lapse(0.002); // [K/m] approximate
45 // Reference: http://en.wikipedia.org/wiki/Lapse_rate
46 return dpsl - dp_lapse * h_ft * atmodel::foot;
49 // Height as a function of pressure.
50 // Valid in the troposphere only.
51 double FGAtmo::a_vs_p(const double press, const double qnh) {
52 using namespace atmodel;
54 double nn = lam0 * Rgas / g / mm;
55 return T0 * ( pow(qnh/P0,nn) - pow(press/P0,nn) ) / lam0;
59 void FGAtmoCache::tabulate() {
60 using namespace atmodel;
62 a_tvs_p = new SGInterpTable;
64 for (double hgt = -1000; hgt <= 32000;) {
65 double press = p_vs_a(hgt);
66 a_tvs_p->addEntry(press / inHg, hgt / foot);
68 #ifdef DEBUG_EXPORT_P_H
70 char* fmt = " { %9.2f , %5.0f },";
71 if (press < 10000) fmt = " { %9.3f , %5.0f },";
72 snprintf(buf, 100, fmt, press, hgt);
83 // make sure cache is valid
84 void FGAtmoCache::cache() {
89 // Pressure within a layer, as a function of height.
90 // Physics model: standard or nonstandard atmosphere,
91 // depending on what parameters you pass in.
92 // Height in meters, pressures in pascals.
93 // As always, lapse is positive in the troposphere,
94 // and zero in the first part of the stratosphere.
96 double FGAtmo::P_layer(const double height, const double href,
97 const double Pref, const double Tref,
99 using namespace atmodel;
101 double N = lapse * Rgas / mm / g;
102 return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
104 return Pref * exp(-g * mm / Rgas / Tref * (height - href));
108 // Check the basic function,
109 // then compare against the interpolator.
110 void FGAtmoCache::check_model() {
129 for (int i = 0; ; i++) {
130 double height = hgts[i];
133 using namespace atmodel;
135 double press = p_vs_a(height);
136 cout << "Height: " << height
137 << " \tpressure: " << press << endl;
139 << a_tvs_p->interpolate(press / inHg)*foot << endl;
143 //////////////////////////////////////////////////////////////////////
145 FGAltimeter::FGAltimeter() : kset(atmodel::ISA::P0), kft(0)
150 double FGAltimeter::reading_ft(const double p_inHg, const double set_inHg) {
151 using namespace atmodel;
152 double press_alt = a_tvs_p->interpolate(p_inHg);
153 double kollsman_shift = a_tvs_p->interpolate(set_inHg);
154 return (press_alt - kollsman_shift);
157 // Altimeter setting.
158 // Field elevation in feet
159 // Field pressure in inHg
160 // field elevation in troposphere only
161 double FGAtmo::qnh(const double field_ft, const double press_inHg) {
162 using namespace atmodel;
164 // Equation derived in altimetry.htm
165 // exponent in QNH equation:
166 double nn = ISA::lam0 * Rgas / g / mm;
167 // pressure ratio factor:
168 double prat = pow(ISA::P0/inHg / press_inHg, nn);
171 * pow(1 + ISA::lam0 * field_ft * foot / ISA::T0 * prat, 1/nn);
174 void FGAltimeter::dump_stack1(const double Tref) {
175 using namespace atmodel;
178 double Psl = P_layer(0, 0, ISA::P0, Tref, ISA::lam0);
179 snprintf(buf, bs, "Tref: %6.2f Psl: %5.0f = %7.4f",
180 Tref, Psl, Psl / inHg);
184 " %6s %6s %6s %6s %6s %6s %6s",
185 "A", "Aind", "Apr", "Aprind", "P", "Psl", "Qnh");
188 double hgts[] = {0, 2500, 5000, 7500, 10000, -9e99};
189 for (int ii = 0; ; ii++) {
190 double hgt_ft = hgts[ii];
193 double press = P_layer(hgt_ft*foot, 0, ISA::P0, Tref, ISA::lam0);
194 double p_inHg = press / inHg;
195 double qnhx = qnh(hgt_ft, p_inHg);
196 double qnh2 = SGMiscd::round(qnhx*100)/100;
198 double Aprind = reading_ft(p_inHg);
199 double Apr = a_vs_p(p_inHg*inHg) / foot;
200 double hind = reading_ft(p_inHg, qnh2);
202 " %6.0f %6.0f %6.0f %6.0f %6.2f %6.2f %6.2f",
203 hgt_ft, hind, Apr, Aprind, p_inHg, Psl/inHg, qnh2);
209 void FGAltimeter::dump_stack() {
210 using namespace atmodel;
211 cout << "........." << endl;
212 cout << "Size: " << sizeof(FGAtmo) << endl;
213 dump_stack1(ISA::T0);
214 dump_stack1(ISA::T0 - 20);