src/Environment/atmosphere.cxx

   1 #include <boost/tuple/tuple.hpp>
   2
   3 #include <simgear/math/SGMath.hxx>
   4 #include <simgear/debug/logstream.hxx>
   5
   6 #include "atmosphere.hxx"
   7
   8 using namespace std;
   9 #include <iostream>
  10 #include <cstdio>
  11
  12 const ISA_layer ISA_def[] = {
  13 //        0    1        2        3           4         5      6         7         8
  14 //       id   (m)      (ft)     (Pa)        (inHg)    (K)     (C)      (K/m)     (K/ft)
  15 ISA_layer(0,      0,       0,   101325,    29.92126, 288.15,  15.00,  0.0065,   0.0019812),
  16 ISA_layer(1,  11000,   36089,  22632.1,    6.683246, 216.65, -56.50,       0,           0),
  17 ISA_layer(2,  20000,   65616,  5474.89,    1.616734, 216.65, -56.50, -0.0010,  -0.0003048),
  18 ISA_layer(3,  32000,  104986,  868.019,    0.256326, 228.65, -44.50, -0.0028,  -0.0008534),
  19 ISA_layer(4,  47000,  154199,  110.906,   0.0327506, 270.65,  -2.50,       0,           0),
  20 ISA_layer(5,  51000,  167322,  66.9389,   0.0197670, 270.65,  -2.50,  0.0028,   0.0008534),
  21 ISA_layer(6,  71000,  232939,  3.95642,  0.00116833, 214.65, -58.50,  0.0020,   0.0006096),
  22 ISA_layer(7,  80000,  262467,  0.88628, 0.000261718, 196.65, -76.50,  0.0,      0.0),
  23 ISA_layer(8,  1.0e9,  3.28e9,  0.00001, 3.0e-9,        2.73, -270.4,  0.0,      0.0),
  24 };
  25
  26 const int ISA_def_size(sizeof(ISA_def) / sizeof(ISA_layer));
  27
  28 // Pressure within a layer, as a function of height.
  29 // Physics model:  standard or nonstandard atmosphere,
  30 //    depending on what parameters you pass in.
  31 // Height in meters, pressures in pascals.
  32 // As always, lapse is positive in the troposphere,
  33 // and zero in the first part of the stratosphere.
  34
  35 double P_layer(const double height, const double href,
  36         const double Pref, const double Tref,
  37         const double lapse) {
  38     using namespace atmodel;
  39     if (lapse) {
  40         double N = lapse * Rgas / mm / g;
  41         return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
  42     } else {
  43         return Pref * exp(-g * mm / Rgas / Tref * (height - href));
  44     }
  45 }
  46
  47
  48 // Temperature within a layer, as a function of height.
  49 // Physics model:  standard or nonstandard atmosphere
  50 //  depending on what parameters you pass in.
  51 // $hh in meters, pressures in Pa.
  52 // As always, $lambda is positive in the troposphere,
  53 // and zero in the first part of the stratosphere.
  54 double T_layer (
  55           const double hh,
  56           const double hb,
  57           const double Pb,
  58           const double Tb,
  59           const double lambda) {
  60   return Tb - lambda*(hh - hb);
  61 }
  62
  63 // Pressure and temperature as a function of height, Psl, and Tsl.
  64 // heights in meters, pressures in Pa.
  65 // Daisy chain version.
  66 // We need "seed" values for sea-level pressure and temperature.
  67 // In addition, for every layer, we need three things
  68 //  from the table: the reference height in that layer,
  69 //  the lapse in that layer, and the cap (if any) for that layer
  70 // (which we take from the /next/ row of the table, if any).
  71 pair<double,double> PT_vs_hpt(
  72       const double hh,
  73       const double _p0,
  74       const double _t0
  75 ) {
  76
  77   const double d0(0);
  78   double hgt = ISA_def[0].height;
  79   double p0 =  _p0;
  80   double t0 =  _t0;
  81 #if 0
  82     cout << "PT_vs_hpt: " << hh << "   " << p0 << "   " << t0 << endl;
  83 #endif
  84
  85   int ii = 0;
  86   for (const ISA_layer* pp = ISA_def; pp->lapse != -1; pp++, ii++) {
  87 #if 0
  88     cout << "PT_vs_hpt: " << ii
  89         << "  height: " << pp->height
  90         << "  temp: "   << pp->temp
  91         << "  lapse: "  << pp->lapse
  92         << endl;
  93 #endif
  94     double xhgt(9e99);
  95     double lapse = pp->lapse;
  96 // Stratosphere starts at a definite temperature,
  97 // not a definite height:
  98     if (ii == 0) {
  99       xhgt = hgt + (t0 - (pp+1)->temp) / lapse;
 100     } else if ((pp+1)->lapse != -1) {
 101       xhgt = (pp+1)->height;
 102     }
 103     if (hh <= xhgt) {
 104       return make_pair(P_layer(hh, hgt, p0, t0, lapse),
 105                  T_layer(hh, hgt, p0, t0, lapse));
 106     }
 107     p0 = P_layer(xhgt, hgt, p0, t0, lapse);
 108     t0 = t0 - lapse * (xhgt - hgt);
 109     hgt = xhgt;
 110   }
 111
 112 // Should never get here.
 113   SG_LOG(SG_ENVIRONMENT, SG_ALERT, "PT_vs_hpt: ran out of layers for h=" << hh );
 114   return make_pair(d0, d0);
 115 }
 116
 117
 118 FGAtmoCache::FGAtmoCache() :
 119     a_tvs_p(0)
 120 {}
 121
 122 FGAtmoCache::~FGAtmoCache() {
 123     delete a_tvs_p;
 124 }
 125
 126
 127 /////////////
 128 // The following two routines are called "fake" because they
 129 // bypass the exceedingly complicated layer model implied by
 130 // the "weather conditioins" popup menu.
 131 // For now we must bypass it for several reasons, including
 132 // the fact that we don't have an "environment" object for
 133 // the airport (only for the airplane).
 134 // degrees C, height in feet
 135 double FGAtmo::fake_T_vs_a_us(const double h_ft,
 136                 const double Tsl) const {
 137     using namespace atmodel;
 138     return Tsl - ISA::lam0 * h_ft * foot;
 139 }
 140
 141 // Dewpoint.  degrees C or K, height in feet
 142 double FGAtmo::fake_dp_vs_a_us(const double dpsl, const double h_ft) {
 143     const double dp_lapse(0.002);       // [K/m] approximate
 144     // Reference:  http://en.wikipedia.org/wiki/Lapse_rate
 145     return dpsl - dp_lapse * h_ft * atmodel::foot;
 146 }
 147
 148 // Height as a function of pressure.
 149 // Valid in the troposphere only.
 150 double FGAtmo::a_vs_p(const double press, const double qnh) {
 151     using namespace atmodel;
 152     using namespace ISA;
 153     double nn = lam0 * Rgas / g / mm;
 154     return T0 * ( pow(qnh/P0,nn) - pow(press/P0,nn) ) / lam0;
 155 }
 156
 157 // force retabulation
 158 void FGAtmoCache::tabulate() {
 159     using namespace atmodel;
 160     delete a_tvs_p;
 161     a_tvs_p = new SGInterpTable;
 162
 163     for (double hgt = -1000; hgt <= 32000;) {
 164         double press,temp;
 165         boost::tie(press, temp) = PT_vs_hpt(hgt);
 166         a_tvs_p->addEntry(press / inHg, hgt / foot);
 167
 168 #ifdef DEBUG_EXPORT_P_H
 169         char buf[100];
 170         char* fmt = " { %9.2f  , %5.0f },";
 171         if (press < 10000) fmt = " {  %9.3f , %5.0f },";
 172         snprintf(buf, 100, fmt, press, hgt);
 173         cout << buf << endl;
 174 #endif
 175         if (hgt < 6000) {
 176             hgt += 500;
 177         } else {
 178             hgt += 1000;
 179         }
 180     }
 181 }
 182
 183 // make sure cache is valid
 184 void FGAtmoCache::cache() {
 185     if (!a_tvs_p)
 186         tabulate();
 187 }
 188
 189 // Check the basic function,
 190 // then compare against the interpolator.
 191 void FGAtmoCache::check_model() {
 192     double hgts[] = {
 193         -1000,
 194         -250,
 195         0,
 196         250,
 197         1000,
 198         5250,
 199         11000,
 200         11000.00001,
 201         15500,
 202         20000,
 203         20000.00001,
 204         25500,
 205         32000,
 206         32000.00001,
 207        -9e99
 208     };
 209
 210     for (int i = 0; ; i++) {
 211         double height = hgts[i];
 212         if (height < -1e6)
 213             break;
 214         using namespace atmodel;
 215         cache();
 216         double press,temp;
 217         boost::tie(press, temp) = PT_vs_hpt(height);
 218         cout << "Height: " << height
 219              << " \tpressure: " << press << endl;
 220         cout << "Check:  "
 221              << a_tvs_p->interpolate(press / inHg)*foot << endl;
 222     }
 223 }
 224
 225 //////////////////////////////////////////////////////////////////////
 226
 227 FGAltimeter::FGAltimeter()
 228 {
 229     cache();
 230 }
 231
 232 double FGAltimeter::reading_ft(const double p_inHg, const double set_inHg) {
 233     using namespace atmodel;
 234     double press_alt      = a_tvs_p->interpolate(p_inHg);
 235     double kollsman_shift = a_tvs_p->interpolate(set_inHg);
 236     return (press_alt - kollsman_shift);
 237 }
 238
 239 // Altimeter setting _in pascals_
 240 //  ... caller gets to convert to inHg or millibars
 241 // Field elevation in m
 242 // Field pressure in pascals
 243 // Valid for fields within the troposphere only.
 244 double FGAtmo::QNH(const double field_elev, const double field_press) {
 245     using namespace atmodel;
 246
 247     //  Equation derived in altimetry.htm
 248     // exponent in QNH equation:
 249     double nn = ISA::lam0 * Rgas / g / mm;
 250     // pressure ratio factor:
 251     double prat = pow(ISA::P0 / field_press, nn);
 252     double rslt =  field_press
 253             * pow(1. + ISA::lam0 * field_elev / ISA::T0 * prat, 1./nn);
 254 #if 0
 255     SG_LOG(SG_ENVIRONMENT, SG_ALERT, "QNH: elev: " << field_elev
 256                 << "  press: " << field_press
 257                 << "  prat: "  << prat
 258                 << "  rslt: " << rslt
 259                 << "  inHg: " << inHg
 260                 << "  rslt/inHG: " << rslt/inHg);
 261 #endif
 262     return rslt;
 263 }
 264
 265 // Invert the QNH calculation to get the field pressure from a metar
 266 // report.
 267 // field pressure _in pascals_
 268 //  ... caller gets to convert to inHg or millibars
 269 // Field elevation in m
 270 // Altimeter setting (QNH) in pascals
 271 // Valid for fields within the troposphere only.
 272 double FGAtmo::fieldPressure(const double field_elev, const double qnh)
 273 {
 274     using namespace atmodel;
 275     static const double nn = ISA::lam0 * Rgas / g / mm;
 276     const double pratio = pow(qnh / ISA::P0, nn);
 277     return ISA::P0 * pow(pratio - field_elev * ISA::lam0 / ISA::T0, 1.0 / nn);
 278 }
 279
 280 void FGAltimeter::dump_stack1(const double Tref) {
 281     using namespace atmodel;
 282     const int bs(200);
 283     char buf[bs];
 284     double Psl = P_layer(0, 0, ISA::P0, Tref, ISA::lam0);
 285     snprintf(buf, bs, "Tref: %6.2f  Psl:  %5.0f = %7.4f",
 286                          Tref,        Psl,     Psl / inHg);
 287     cout << buf << endl;
 288
 289     snprintf(buf, bs,
 290         " %6s  %6s  %6s  %6s   %6s  %6s   %6s",
 291         "A", "Aind", "Apr", "Aprind", "P", "Psl", "Qnh");
 292     cout << buf << endl;
 293
 294     double hgts[] = {0, 2500, 5000, 7500, 10000, -9e99};
 295     for (int ii = 0; ; ii++) {
 296         double hgt_ft = hgts[ii];
 297         double hgt = hgt_ft * foot;
 298         if (hgt_ft < -1e6)
 299             break;
 300         double press = P_layer(hgt, 0, ISA::P0, Tref, ISA::lam0);
 301         double qnhx = QNH(hgt, press) / inHg;
 302         double qnh2 = SGMiscd::round(qnhx*100)/100;
 303
 304         double p_inHg = press / inHg;
 305         double Aprind = reading_ft(p_inHg);
 306         double Apr    = a_vs_p(p_inHg*inHg) / foot;
 307         double hind = reading_ft(p_inHg, qnh2);
 308         snprintf(buf, bs,
 309             " %6.0f  %6.0f  %6.0f  %6.0f   %6.2f  %6.2f   %6.2f",
 310             hgt_ft,  hind,   Apr,  Aprind,  p_inHg, Psl/inHg, qnh2);
 311         cout << buf << endl;
 312     }
 313 }
 314
 315
 316 void FGAltimeter::dump_stack() {
 317     using namespace atmodel;
 318     cout << "........." << endl;
 319     cout << "Size: " << sizeof(FGAtmo) << endl;
 320     dump_stack1(ISA::T0);
 321     dump_stack1(ISA::T0 - 20);
 322 }