src/Environment/atmosphere.cxx

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