1 #include "atmosphere.hxx"
6 FGAtmoCache::FGAtmoCache() :
10 FGAtmoCache::~FGAtmoCache() {
14 // Pressure as a function of height.
15 // Valid below 32000 meters,
16 // i.e. troposphere and first two layers of stratosphere.
17 // Does not depend on any caching; can be used to
18 // *construct* caches and interpolation tables.
20 // Height in meters, pressure in pascals.
22 double FGAtmo::p_vs_a(const double height) {
23 using namespace atmodel;
24 if (height <= 11000.) {
25 return P_layer(height, 0.0, ISA::P0, ISA::T0, ISA::lam0);
26 } else if (height <= 20000.) {
27 return P_layer(height, 11000., 22632.06, 216.65, 0.0);
28 } else if (height <= 32000.) {
29 return P_layer(height, 20000., 5474.89, 216.65, -0.001);
34 // degrees C, height in feet
35 double FGAtmo::fake_t_vs_a_us(const double h_ft) {
36 using namespace atmodel;
37 return ISA::T0 - ISA::lam0 * h_ft * foot - freezing;
40 // Dewpoint. degrees C or K, height in feet
41 double FGAtmo::fake_dp_vs_a_us(const double dpsl, const double h_ft) {
42 const double dp_lapse(0.002); // [K/m] approximate
43 // Reference: http://en.wikipedia.org/wiki/Lapse_rate
44 return dpsl - dp_lapse * h_ft * atmodel::foot;
47 // Height as a function of pressure.
48 // Valid in the troposphere only.
49 double FGAtmo::a_vs_p(const double press, const double qnh) {
50 using namespace atmodel;
52 double nn = lam0 * Rgas / g / mm;
53 return T0 * ( pow(qnh/P0,nn) - pow(press/P0,nn) ) / lam0;
57 void FGAtmoCache::tabulate() {
58 using namespace atmodel;
60 a_tvs_p = new SGInterpTable;
62 for (double hgt = -1000; hgt <= 32000;) {
63 double press = p_vs_a(hgt);
64 a_tvs_p->addEntry(press / inHg, hgt / foot);
66 #ifdef DEBUG_EXPORT_P_H
68 char* fmt = " { %9.2f , %5.0f },";
69 if (press < 10000) fmt = " { %9.3f , %5.0f },";
70 snprintf(buf, 100, fmt, press, hgt);
81 // make sure cache is valid
82 void FGAtmoCache::cache() {
87 // Pressure within a layer, as a function of height.
88 // Physics model: standard or nonstandard atmosphere,
89 // depending on what parameters you pass in.
90 // Height in meters, pressures in pascals.
91 // As always, lapse is positive in the troposphere,
92 // and zero in the first part of the stratosphere.
94 double FGAtmo::P_layer(const double height, const double href,
95 const double Pref, const double Tref,
97 using namespace atmodel;
99 double N = lapse * Rgas / mm / g;
100 return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
102 return Pref * exp(-g * mm / Rgas / Tref * (height - href));
106 // Check the basic function,
107 // then compare against the interpolator.
108 void FGAtmoCache::check_model() {
127 for (int i = 0; ; i++) {
128 double height = hgts[i];
131 using namespace atmodel;
133 double press = p_vs_a(height);
134 cout << "Height: " << height
135 << " \tpressure: " << press << endl;
137 << a_tvs_p->interpolate(press / inHg)*foot << endl;
141 //////////////////////////////////////////////////////////////////////
143 FGAltimeter::FGAltimeter() : kset(atmodel::ISA::P0), kft(0)
148 double FGAltimeter::reading_ft(const double p_inHg, const double set_inHg) {
149 using namespace atmodel;
150 double press_alt = a_tvs_p->interpolate(p_inHg);
151 double kollsman_shift = a_tvs_p->interpolate(set_inHg);
152 return (press_alt - kollsman_shift);
155 // Altimeter setting.
156 // Field elevation in feet
157 // Field pressure in inHg
158 // field elevation in troposphere only
159 double FGAtmo::qnh(const double field_ft, const double press_inHg) {
160 using namespace atmodel;
162 // Equation derived in altimetry.htm
163 // exponent in QNH equation:
164 double nn = ISA::lam0 * Rgas / g / mm;
165 // pressure ratio factor:
166 double prat = pow(ISA::P0/inHg / press_inHg, nn);
169 * pow(1 + ISA::lam0 * field_ft * foot / ISA::T0 * prat, 1/nn);
172 void FGAltimeter::dump_stack1(const double Tref) {
173 using namespace atmodel;
176 double Psl = P_layer(0, 0, ISA::P0, Tref, ISA::lam0);
177 snprintf(buf, bs, "Tref: %6.2f Psl: %5.0f = %7.4f",
178 Tref, Psl, Psl / inHg);
182 " %6s %6s %6s %6s %6s %6s %6s",
183 "A", "Aind", "Apr", "Aprind", "P", "Psl", "Qnh");
186 double hgts[] = {0, 2500, 5000, 7500, 10000, -9e99};
187 for (int ii = 0; ; ii++) {
188 double hgt_ft = hgts[ii];
191 double press = P_layer(hgt_ft*foot, 0, ISA::P0, Tref, ISA::lam0);
192 double p_inHg = press / inHg;
193 double qnhx = qnh(hgt_ft, p_inHg);
194 double qnh2 = round(qnhx*100)/100;
196 double Aprind = reading_ft(p_inHg);
197 double Apr = a_vs_p(p_inHg*inHg) / foot;
198 double hind = reading_ft(p_inHg, qnh2);
200 " %6.0f %6.0f %6.0f %6.0f %6.2f %6.2f %6.2f",
201 hgt_ft, hind, Apr, Aprind, p_inHg, Psl/inHg, qnh2);
207 void FGAltimeter::dump_stack() {
208 using namespace atmodel;
209 cout << "........." << endl;
210 cout << "Size: " << sizeof(FGAtmo) << endl;
211 dump_stack1(ISA::T0);
212 dump_stack1(ISA::T0 - 20);