+const ISA_layer ISA_def[] = {
+// 0 1 2 3 4 5 6 7 8
+// id (m) (ft) (Pa) (inHg) (K) (C) (K/m) (K/ft)
+ISA_layer(0, 0, 0, 101325, 29.92126, 288.15, 15.00, 0.0065, 0.0019812),
+ISA_layer(1, 11000, 36089, 22632.1, 6.683246, 216.65, -56.50, 0, 0),
+ISA_layer(2, 20000, 65616, 5474.89, 1.616734, 216.65, -56.50, -0.0010, -0.0003048),
+ISA_layer(3, 32000, 104986, 868.019, 0.256326, 228.65, -44.50, -0.0028, -0.0008534),
+ISA_layer(4, 47000, 154199, 110.906, 0.0327506, 270.65, -2.50, 0, 0),
+ISA_layer(5, 51000, 167322, 66.9389, 0.0197670, 270.65, -2.50, 0.0028, 0.0008534),
+ISA_layer(6, 71000, 232939, 3.95642, 0.00116833, 214.65, -58.50, 0.0020, 0.0006096),
+ISA_layer(7, 80000, 262467, 0.88628, 0.000261718, 196.65, -76.50),
+};
+
+const int ISA_def_size(sizeof(ISA_def) / sizeof(ISA_layer));
+
+// Pressure within a layer, as a function of height.
+// Physics model: standard or nonstandard atmosphere,
+// depending on what parameters you pass in.
+// Height in meters, pressures in pascals.
+// As always, lapse is positive in the troposphere,
+// and zero in the first part of the stratosphere.
+
+double P_layer(const double height, const double href,
+ const double Pref, const double Tref,
+ const double lapse) {
+ using namespace atmodel;
+ if (lapse) {
+ double N = lapse * Rgas / mm / g;
+ return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
+ } else {
+ return Pref * exp(-g * mm / Rgas / Tref * (height - href));
+ }
+}
+
+
+// Temperature within a layer, as a function of height.
+// Physics model: standard or nonstandard atmosphere
+// depending on what parameters you pass in.
+// $hh in meters, pressures in Pa.
+// As always, $lambda is positive in the troposphere,
+// and zero in the first part of the stratosphere.
+double T_layer (
+ const double hh,
+ const double hb,
+ const double Pb,
+ const double Tb,
+ const double lambda) {
+ return Tb - lambda*(hh - hb);
+}
+
+// Pressure and temperature as a function of height, Psl, and Tsl.
+// heights in meters, pressures in Pa.
+// Daisy chain version.
+// We need "seed" values for sea-level pressure and temperature.
+// In addition, for every layer, we need three things
+// from the table: the reference height in that layer,
+// the lapse in that layer, and the cap (if any) for that layer
+// (which we take from the /next/ row of the table, if any).
+pair<double,double> PT_vs_hpt(
+ const double hh,
+ const double _p0,
+ const double _t0
+) {
+
+ const double d0(0);
+ double hgt = ISA_def[0].height;
+ double p0 = _p0;
+ double t0 = _t0;
+#if 0
+ cout << "PT_vs_hpt: " << hh << " " << p0 << " " << t0 << endl;
+#endif
+
+ int ii = 0;
+ for (const ISA_layer* pp = ISA_def; pp->lapse != -1; pp++, ii++) {
+#if 0
+ cout << "PT_vs_hpt: " << ii
+ << " height: " << pp->height
+ << " temp: " << pp->temp
+ << " lapse: " << pp->lapse
+ << endl;
+#endif
+ double xhgt(9e99);
+ double lapse = pp->lapse;
+// Stratosphere starts at a definite temperature,
+// not a definite height:
+ if (ii == 0) {
+ xhgt = hgt + (t0 - (pp+1)->temp) / lapse;
+ } else if ((pp+1)->lapse != -1) {
+ xhgt = (pp+1)->height;
+ }
+ if (hh <= xhgt) {
+ return make_pair(P_layer(hh, hgt, p0, t0, lapse),
+ T_layer(hh, hgt, p0, t0, lapse));
+ }
+ p0 = P_layer(xhgt, hgt, p0, t0, lapse);
+ t0 = t0 - lapse * (xhgt - hgt);
+ hgt = xhgt;
+ }
+
+// Should never get here.
+ SG_LOG(SG_GENERAL, SG_ALERT, "PT_vs_hpt: ran out of layers");
+ return make_pair(d0, d0);
+}
+
+