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