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