src/Environment/atmosphere.cxx

   1 #include <simgear/math/SGMath.hxx>
   2
   3 #include "atmosphere.hxx"
   4
   5 using namespace std;
   6 #include <iostream>
   7
   8 FGAtmoCache::FGAtmoCache() :
   9     a_tvs_p(0)
  10 {}
  11
  12 FGAtmoCache::~FGAtmoCache() {
  13     delete a_tvs_p;
  14 }
  15
  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.
  21 //
  22 // Height in meters, pressure in pascals.
  23
  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);
  32     }
  33     return 0;
  34 }
  35
  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;
  40 }
  41
  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;
  47 }
  48
  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;
  53     using namespace ISA;
  54     double nn = lam0 * Rgas / g / mm;
  55     return T0 * ( pow(qnh/P0,nn) - pow(press/P0,nn) ) / lam0;
  56 }
  57
  58 // force retabulation
  59 void FGAtmoCache::tabulate() {
  60     using namespace atmodel;
  61     delete a_tvs_p;
  62     a_tvs_p = new SGInterpTable;
  63
  64     for (double hgt = -1000; hgt <= 32000;) {
  65         double press = p_vs_a(hgt);
  66         a_tvs_p->addEntry(press / inHg, hgt / foot);
  67
  68 #ifdef DEBUG_EXPORT_P_H
  69         char buf[100];
  70         char* fmt = " { %9.2f  , %5.0f },";
  71         if (press < 10000) fmt = " {  %9.3f , %5.0f },";
  72         snprintf(buf, 100, fmt, press, hgt);
  73         cout << buf << endl;
  74 #endif
  75         if (hgt < 6000) {
  76             hgt += 500;
  77         } else {
  78             hgt += 1000;
  79         }
  80     }
  81 }
  82
  83 // make sure cache is valid
  84 void FGAtmoCache::cache() {
  85     if (!a_tvs_p)
  86         tabulate();
  87 }
  88
  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.
  95
  96 double FGAtmo::P_layer(const double height, const double href,
  97         const double Pref, const double Tref,
  98         const double lapse) {
  99     using namespace atmodel;
 100     if (lapse) {
 101         double N = lapse * Rgas / mm / g;
 102         return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
 103     } else {
 104         return Pref * exp(-g * mm / Rgas / Tref * (height - href));
 105     }
 106 }
 107
 108 // Check the basic function,
 109 // then compare against the interpolator.
 110 void FGAtmoCache::check_model() {
 111     double hgts[] = {
 112         -1000,
 113         -250,
 114         0,
 115         250,
 116         1000,
 117         5250,
 118         11000,
 119         11000.00001,
 120         15500,
 121         20000,
 122         20000.00001,
 123         25500,
 124         32000,
 125         32000.00001,
 126        -9e99
 127     };
 128
 129     for (int i = 0; ; i++) {
 130         double height = hgts[i];
 131         if (height < -1e6)
 132             break;
 133         using namespace atmodel;
 134         cache();
 135         double press = p_vs_a(height);
 136         cout << "Height: " << height
 137              << " \tpressure: " << press << endl;
 138         cout << "Check:  "
 139              << a_tvs_p->interpolate(press / inHg)*foot << endl;
 140     }
 141 }
 142
 143 //////////////////////////////////////////////////////////////////////
 144
 145 FGAltimeter::FGAltimeter() : kset(atmodel::ISA::P0), kft(0)
 146 {
 147     cache();
 148 }
 149
 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);
 155 }
 156
 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;
 163
 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);
 169
 170     return press_inHg
 171             * pow(1 + ISA::lam0 * field_ft * foot / ISA::T0 * prat, 1/nn);
 172 }
 173
 174 void FGAltimeter::dump_stack1(const double Tref) {
 175     using namespace atmodel;
 176     const int bs(200);
 177     char buf[bs];
 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);
 181     cout << buf << endl;
 182
 183     snprintf(buf, bs,
 184         " %6s  %6s  %6s  %6s   %6s  %6s   %6s",
 185         "A", "Aind", "Apr", "Aprind", "P", "Psl", "Qnh");
 186     cout << buf << endl;
 187
 188     double hgts[] = {0, 2500, 5000, 7500, 10000, -9e99};
 189     for (int ii = 0; ; ii++) {
 190         double hgt_ft = hgts[ii];
 191         if (hgt_ft < -1e6)
 192             break;
 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;
 197
 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);
 201         snprintf(buf, bs,
 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);
 204         cout << buf << endl;
 205     }
 206 }
 207
 208
 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);
 215 }