src/Environment/atmosphere.cxx

   1 #include "atmosphere.hxx"
   2
   3 using namespace std;
   4 #include <iostream>
   5
   6 FGAtmoCache::FGAtmoCache() :
   7     a_tvs_p(0)
   8 {}
   9
  10 FGAtmoCache::~FGAtmoCache() {
  11     delete a_tvs_p;
  12 }
  13
  14 // Pressure as a function of height.
  15 // Valid below 32000 meters,
  16 //   i.e. troposphere and first two layers of stratosphere.
  17 // Does not depend on any caching; can be used to
  18 // *construct* caches and interpolation tables.
  19 //
  20 // Height in meters, pressure in pascals.
  21
  22 double FGAtmo::p_vs_a(const double height) {
  23     using namespace atmodel;
  24     if (height <= 11000.) {
  25         return P_layer(height, 0.0, ISA::P0, ISA::T0, ISA::lam0);
  26     } else if (height <= 20000.) {
  27         return P_layer(height, 11000., 22632.06, 216.65, 0.0);
  28     } else if (height <= 32000.) {
  29         return P_layer(height, 20000.,  5474.89, 216.65, -0.001);
  30     }
  31     return 0;
  32 }
  33
  34 // degrees C, height in feet
  35 double FGAtmo::fake_t_vs_a_us(const double h_ft) {
  36     using namespace atmodel;
  37     return ISA::T0 - ISA::lam0 * h_ft * foot - freezing;
  38 }
  39
  40 // Dewpoint.  degrees C or K, height in feet
  41 double FGAtmo::fake_dp_vs_a_us(const double dpsl, const double h_ft) {
  42     const double dp_lapse(0.002);       // [K/m] approximate
  43     // Reference:  http://en.wikipedia.org/wiki/Lapse_rate
  44     return dpsl - dp_lapse * h_ft * atmodel::foot;
  45 }
  46
  47 // Height as a function of pressure.
  48 // Valid in the troposphere only.
  49 double FGAtmo::a_vs_p(const double press, const double qnh) {
  50     using namespace atmodel;
  51     using namespace ISA;
  52     double nn = lam0 * Rgas / g / mm;
  53     return T0 * ( pow(qnh/P0,nn) - pow(press/P0,nn) ) / lam0;
  54 }
  55
  56 // force retabulation
  57 void FGAtmoCache::tabulate() {
  58     using namespace atmodel;
  59     delete a_tvs_p;
  60     a_tvs_p = new SGInterpTable;
  61
  62     for (double hgt = -1000; hgt <= 32000;) {
  63         double press = p_vs_a(hgt);
  64         a_tvs_p->addEntry(press / inHg, hgt / foot);
  65
  66 #ifdef DEBUG_EXPORT_P_H
  67         char buf[100];
  68         char* fmt = " { %9.2f  , %5.0f },";
  69         if (press < 10000) fmt = " {  %9.3f , %5.0f },";
  70         snprintf(buf, 100, fmt, press, hgt);
  71         cout << buf << endl;
  72 #endif
  73         if (hgt < 6000) {
  74             hgt += 500;
  75         } else {
  76             hgt += 1000;
  77         }
  78     }
  79 }
  80
  81 // make sure cache is valid
  82 void FGAtmoCache::cache() {
  83     if (!a_tvs_p)
  84         tabulate();
  85 }
  86
  87 // Pressure within a layer, as a function of height.
  88 // Physics model:  standard or nonstandard atmosphere,
  89 //    depending on what parameters you pass in.
  90 // Height in meters, pressures in pascals.
  91 // As always, lapse is positive in the troposphere,
  92 // and zero in the first part of the stratosphere.
  93
  94 double FGAtmo::P_layer(const double height, const double href,
  95         const double Pref, const double Tref,
  96         const double lapse) {
  97     using namespace atmodel;
  98     if (lapse) {
  99         double N = lapse * Rgas / mm / g;
 100         return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
 101     } else {
 102         return Pref * exp(-g * mm / Rgas / Tref * (height - href));
 103     }
 104 }
 105
 106 // Check the basic function,
 107 // then compare against the interpolator.
 108 void FGAtmoCache::check_model() {
 109     double hgts[] = {
 110         -1000,
 111         -250,
 112         0,
 113         250,
 114         1000,
 115         5250,
 116         11000,
 117         11000.00001,
 118         15500,
 119         20000,
 120         20000.00001,
 121         25500,
 122         32000,
 123         32000.00001,
 124        -9e99
 125     };
 126
 127     for (int i = 0; ; i++) {
 128         double height = hgts[i];
 129         if (height < -1e6)
 130             break;
 131         using namespace atmodel;
 132         cache();
 133         double press = p_vs_a(height);
 134         cout << "Height: " << height
 135              << " \tpressure: " << press << endl;
 136         cout << "Check:  "
 137              << a_tvs_p->interpolate(press / inHg)*foot << endl;
 138     }
 139 }
 140
 141 //////////////////////////////////////////////////////////////////////
 142
 143 FGAltimeter::FGAltimeter() : kset(atmodel::ISA::P0), kft(0)
 144 {
 145     cache();
 146 }
 147
 148 double FGAltimeter::reading_ft(const double p_inHg, const double set_inHg) {
 149     using namespace atmodel;
 150     double press_alt      = a_tvs_p->interpolate(p_inHg);
 151     double kollsman_shift = a_tvs_p->interpolate(set_inHg);
 152     return (press_alt - kollsman_shift);
 153 }
 154
 155 // Altimeter setting.
 156 // Field elevation in feet
 157 // Field pressure in inHg
 158 // field elevation in troposphere only
 159 double FGAtmo::qnh(const double field_ft, const double press_inHg) {
 160     using namespace atmodel;
 161
 162     //  Equation derived in altimetry.htm
 163     // exponent in QNH equation:
 164     double nn = ISA::lam0 * Rgas / g / mm;
 165     // pressure ratio factor:
 166     double prat = pow(ISA::P0/inHg / press_inHg, nn);
 167
 168     return press_inHg
 169             * pow(1 + ISA::lam0 * field_ft * foot / ISA::T0 * prat, 1/nn);
 170 }
 171
 172 void FGAltimeter::dump_stack1(const double Tref) {
 173     using namespace atmodel;
 174     const int bs(200);
 175     char buf[bs];
 176     double Psl = P_layer(0, 0, ISA::P0, Tref, ISA::lam0);
 177     snprintf(buf, bs, "Tref: %6.2f  Psl:  %5.0f = %7.4f",
 178                          Tref,        Psl,     Psl / inHg);
 179     cout << buf << endl;
 180
 181     snprintf(buf, bs,
 182         " %6s  %6s  %6s  %6s   %6s  %6s   %6s",
 183         "A", "Aind", "Apr", "Aprind", "P", "Psl", "Qnh");
 184     cout << buf << endl;
 185
 186     double hgts[] = {0, 2500, 5000, 7500, 10000, -9e99};
 187     for (int ii = 0; ; ii++) {
 188         double hgt_ft = hgts[ii];
 189         if (hgt_ft < -1e6)
 190             break;
 191         double press = P_layer(hgt_ft*foot, 0, ISA::P0, Tref, ISA::lam0);
 192         double p_inHg = press / inHg;
 193         double qnhx = qnh(hgt_ft, p_inHg);
 194         double qnh2 = round(qnhx*100)/100;
 195
 196         double Aprind = reading_ft(p_inHg);
 197         double Apr    = a_vs_p(p_inHg*inHg) / foot;
 198         double hind = reading_ft(p_inHg, qnh2);
 199         snprintf(buf, bs,
 200             " %6.0f  %6.0f  %6.0f  %6.0f   %6.2f  %6.2f   %6.2f",
 201             hgt_ft,  hind,   Apr,  Aprind,  p_inHg, Psl/inHg, qnh2);
 202         cout << buf << endl;
 203     }
 204 }
 205
 206
 207 void FGAltimeter::dump_stack() {
 208     using namespace atmodel;
 209     cout << "........." << endl;
 210     cout << "Size: " << sizeof(FGAtmo) << endl;
 211     dump_stack1(ISA::T0);
 212     dump_stack1(ISA::T0 - 20);
 213 }