#! /usr/bin/perl -CS
## Calculate local barometric pressure, based on nearby weather reports.

use strict;
use warnings;
use ArgParser;
use utf8;

my $usage_msg =  < it strictly follows the FAA algorithm for AWOS.
> The algorithm for the Altimeter is attached for your information.

I quote from the attachment:  In normal operation:

   AS = (Pa^N + K Ha)^(1/N)
where:
      Pa = Field pressure
      N = 0.1903
      K = 1.313e-5
      Ha = Field elevation of the airport
           in feet above mean sea level

This algorithm creates a one-to-one mapping between
field pressure and the reported altimeter setting.
If all we want is field pressure, that is all we need to
know.

However, you can make a plausibility argument in favor of
those numbers by showing that in the International Standard
Atmosphere, different stations will agree as to the
altimeter setting (to a fairly good approximation),
independent of the altitude of the station.  This gives
altimeter settings fairly good "portability" -- which is
highly desirable in the aeronautical application for which
altimeter settings are intended.

For more discussion of the above formula, and its connection
to altimetry and to the standard atmosphere, see
  http://www.av8n.com/physics/altimetry.htm
EoF

my $baro_details = <new(\@ARGV);
  sub check {return $scan->check(@_)};

  set_nature();

  while (@ARGV) {
    my $verb = $scan->get_arg();
    $scan->{arg} =~ s'^--'-';
    if (0) {}
    elsif (check('-hh|-H')) {
      print $baro_msg;
      print $baro_details;
      exit(0)
    }

    elsif (check('-h')) {
      print $usage_msg;
      exit(0);
    }

    elsif (check('-v')) {
      $X::verbose++;
    }
    elsif (check('-formula')) {
      print $formula;
    }
    elsif (check('-scale')) {
      scale();

    }
    elsif (check('-halt|-HALT')) {
      HALT();
    }
    elsif (check('-check1', 'p__inHg, T__C')) {
      my $arg0 = $scan->get_word() * $X::inHg;
      my $arg1 = $scan->get_word() + $X::freezing;
      check_PhTh($arg0, $arg1);
    }
    elsif (check('-tabulate' )) {
      tabulate_Ph();
    }
    elsif (check('-tex' )) {
      tex();
    }
    elsif (check('-cxx' )) {
      cxx();
    }
    elsif (check('-check2' )) {
      check2();
    }
    elsif (check('-sweep', 'h__feet')) {
      my $hmax = $scan->get_word();
      do_sweep($hmax);
    }
    elsif (check('-baro')) {
      baro();
    }
    elsif (check('-interp')) {
      interpolate_Ph_fgfs();         ## for FlightGear
    }
    elsif (check '-show-cmds') {
      $scan->show_cmds;
    }
    elsif ($scan->end_learning) {
    }
    else {
      die "Extraneous verbiage: '$verb'\n";
    }
  }
  exit(0);
}

## Pressure as a function of height in the standard atmosphere.
## In other words, pressure as a function of pressure_altitude
## heights in meters, pressures in Pa.
## Daisy chain version.
## For the first layer, we need initial values for
##  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 for that layer (which we take from
##  the /next/ row of the table, if any).
sub PhTh {
  my ($hh, $psl, $tsl) = @_;
  my $maxrow = $#X::ISA_table;
  my $row = $X::ISA_table[0];
  my ($ndx, $hgt, $hgt_ft, $p0, $p0_in, $t0, $lapse, $lapse_ft) = @$row;
  if (defined $psl) {
    $p0 = $psl;
  }
  if (defined $tsl) {
    $t0 = $tsl;
  }
  for (my $ii = 0; $ii <= $maxrow; $ii++){
    my $xhgt;
    my $xrow;
    $xrow = $X::ISA_table[1+$ii];
## Stratosphere starts at a definite temperature,
## not a definite height:
    if ($ii == 0) {
      $xhgt = $hgt + ($t0 - $xrow->[5]) / $lapse;
    }
    elsif ($ii < $maxrow) {
      $xhgt = $xrow->[1];
    }
    if (!defined $xhgt || $hh <= $xhgt ) {
      return (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;
    $lapse = $xrow->[6];
  }

  die "Ran out of layers";    # should never happen
}


sub cxx {
#           0    1        2        3           4         5      6         7         8
#          id   (m)      (ft)     (Pa)        (inHg)    (K)     (C)      (K/m)     (K/ft)
# ISA_layer(2,  20000,   65616,  5474.89,    1.616734, 216.65, -56.50, -0.0010,  -0.0003048),
  my @todo = (0, 11000, 20000, 32000, 47000, 51000, 71000, 80000);

  my $jj = 0;
  foreach my $alt (@todo) {
    my ($ii, $altx, $alt_ftx, $press, $press_inx, $tempK,  $tempC,
          $lapsex, $lapse_ftx) = fmt_row($alt, $jj);

    $altx =~ s','';
    $alt_ftx =~ s','';
    printf('ISA_layer(%d, %6s, %7s, %8.6g, %11s, %6.2f, %6.2f',
        $ii, $altx, $alt_ftx, $press, $press_inx, $tempK,  $tempC);

    if (defined $lapsex) {
      $lapsex =~ s'~'';
      $lapse_ftx =~ s'~'';
      printf(', %7s,  %10s', $lapsex, $lapse_ftx);
    }

    print "),\n";
    $jj++;
  }
}

sub tex {
#         (m)      (ft)       (Pa)          (inHg)      m/kg^3  (K)       (C)     boundary     num      (K/m)     (K/ft)
# \X{}{ 11,000 }{  36,089 }{  22632.1 }{     6.68325 }{ xxx }{ 216.65 }{ -56.50 }{Tropopause}\Y{1}{      ~0 }{         ~0 }{Stratosphere}

  my @todo     = (0, 11000, 20000, 32000,    47000,   51000, 71000, 80000);
  my @boundary = ("Sea Level", "Tropopause", "",  "", "Stratopause", "", "", "Mesopause");
  my @interior = ("Troposphere", "Stratosphere","Stratosphere","Stratosphere",
                  "Mesosphere","Mesosphere","Mesosphere");

  my $jj = 0;
  foreach my $alt (@todo) {
    my ($ii, $altx, $alt_ftx, $press, $press_inx, $tempK,  $tempC,
          $lapsex, $lapse_ftx) = fmt_row($alt, $jj);

# pv = n kt  n/v = p/kt

    my $density = $X::mm*$press/$tempK/$X::Rgas;

    printf('\X{}{ %6s }{ %7s }{ %8.6g }{ %11s }{ %8.6f }{ %6.2f }{ %6.2f }{ %11s }',
           $altx, $alt_ftx, $press, $press_inx, $density, $tempK,  $tempC, $boundary[$ii]);

    if (defined $lapsex) {
      printf('\Y{%d}{ %7s }{ %10s }{ %11s }', $ii, $lapsex, $lapse_ftx, $interior[$ii]);
    }

    print "\n";
    $jj++;
  }
}

sub fmt_row {
  my ($alt, $ii) = @_;
  my @rslt;

    my $lapse = ${$X::ISA_table[$ii]}[6];
    my ($press, $temp) = PhTh($alt);
    my $altx = sprintf("%d", $alt);
    $altx =~ s'000$',000';
    my $alt_ftx = sprintf("%d", $alt / $X::foot);
    $alt_ftx =~ s/(...)$/,$1/;
    my $press_in = $press/$X::inHg;
    my $press_inx = sprintf("%11.7g", $press_in);
    if ($press_in < 10) {
      my $ndig = 6 - int(log($press_in)/log(10));
      $press_inx = sprintf("%11." . $ndig . "f", $press_in);
    }
    if (0 && $press_in < 1) {
      $press_inx = sprintf("%11.5f", $press_in);
    }
    @rslt = ($ii, $altx, $alt_ftx, $press, $press_inx,
          $temp,$temp-$X::freezing);

    if (defined $lapse) {
      my $lapsex = sprintf("%+7.4f", $lapse);
      my $lapse_ftx = sprintf("%+10.7f", $lapse*$X::foot);
      $lapsex =~ s'0[.]0*$'0';
      $lapsex =~ s'[+]'~';
      $lapse_ftx =~ s'0[.]0*$'0';
      $lapse_ftx =~ s'[+]'~';
      push @rslt, ($lapsex, $lapse_ftx);
    }
    return @rslt;
}


sub check2 {
# 0    1      2          3          4             5        6            7
#(b)  (m)    (ft)       (pascals)  (inHg)        (K)      (K/m)        (K/ft)
  my $row0 = $X::ISA_table[0];
  if (@X::ISA_table_wiki_orig) {}
  my $t0 = $row0->[5] - $X::freezing;
  ($t0 == 15.) || die "bad t0 in table\n";
  ################     b   m   fg   Pa     inHg
  my $fmt  = "table:   %d %5d  %7d %9.6g %12.8f %8.2f %8.4f %12.8f\n";
  my $xfmt = "predict: %d %5d  %7d %9.6g %12.8f %8.2f\n";

#  foreach my $row (@X::ISA_table) {
#    printf($fmt, @$row);
#  }

  my $maxrow = $#X::ISA_table;
  for (my $ii = 0; $ii <= $maxrow; $ii++){
    my $row = $X::ISA_table[$ii];
    my ($ndx, $hgt, $hgt_ft, $p0, $p0_in, $t0, $lapse, $lapse_ft) = @$row;
    printf($fmt, $ndx, $hgt, $hgt_ft, $p0, $p0_in, $t0, $lapse, $lapse_ft);

    my ($xndx, $xhgt, $xhgt_ft, $xp0, $xp0_in, $xt0, $xlapse, $xlapse_ft);
    print "\n";
    if ($ii < $maxrow) {
      my $xrow = $X::ISA_table[1+$ii];
      ($xndx, $xhgt, $xhgt_ft, $xp0, $xp0_in, $xt0, $xlapse, $xlapse_ft) = @$xrow;
    } else {
      $xndx = 1+$maxrow;
      $xhgt = 80000;
      $xhgt_ft = $xhgt / $X::foot;
    }
#     $xp0 =  P_layer($xhgt, $hgt, $p0, $t0, $lapse);
    $xp0 =  P_layer($xhgt, $hgt, $p0_in*$X::inHg, $t0, $lapse);
    $xp0_in = $xp0 / $X::inHg;
    $xt0 =  T_layer($xhgt, $hgt, $p0, $t0, $lapse);
    printf($xfmt, $xndx, $xhgt, $xhgt_ft, $xp0, $xp0_in, $xt0);
  }
}

# Check pressure and temperature.
# Check for continuity across layer boundaries:
sub check_PhTh {
  my ($psl, $tsl) = @_;
  my @todo = (-1000, 0, 1000,
                      11000, 11000.00001,
                        12000,
                      20000, 20000.00001,
                      32000, 32000.00001,
                      47000, 47000.00001,
                      51000, 51000.00001,
                      71000, 71000.00001,
                      80000
                      );
  my $previous = -9e99;
  for my $hgt (sort @todo) {
    my $hgtx = sprintf("%12.5f", $hgt);
    print "hgt: $hgtx   ";
    my ($press,$temp) = PhTh($hgt, $psl, $tsl);
    $press = sprintf("%9.6g", $press);
    $temp = sprintf("%9.3f", $temp);
    if ($press > -10000) {
      print  "press: $press   ";
      $press = sprintf("%9.2f", $press/$X::inHg);
      $hgtx = sprintf("%6d", $hgt/$X::foot);
      print "hgt/ft: $hgtx  press/inHg: $press";
      print " temp: $temp"
    }
    print "\n";
    if (abs($hgt - $previous) < 1) {
      print "\n";
    }
    $previous = $hgt;
  }
  my ($p80,$t80) = PhTh(80000, $psl, $tsl);
  my $pref80 = 0.886279504097685;
  my $delta = ($p80 - $pref80) / $pref80;
  printf("Regression: %8.2g\n", $delta);
}



sub set_nature {
@X::ISA_table = (
#vvv  ............... boundary ......................... vvvv layer above bdy vvvvvvv
# 0    1      2          3                4             5        6            7
#(b)  (m)    (ft)       (pascals)        (inHg)        (K)      (K/m)        (K/ft)
[0,   0,       0,         101325,        29.92126,     288.15,  0.0065,      0.0019812   ],
[1,   11000,  36089,      22632.06409,   6.683245,     216.65,  0.0,         0.0         ],
[2,   20000,  65617,      5474.88813,    1.616734,     216.65,  -0.001,      -0.0003048  ],
[3,   32000,  104987,     868.01872,     0.2563258,    228.65,  -0.0028,     -0.00085344 ],
[4,   47000,  154199,     110.90630,     0.0327506,    270.65,  0.0,         0.0         ],
[5,   51000,  167323,     66.93884,      0.01976704,   270.65,  0.0028,      0.00085344  ],
[6,   71000,  232940,     3.95642,       0.00116833,   214.65,  0.002,       0.0006097   ],
);

## Some constants of nature,
##  and standard units of measure:
  $X::g = 9.80665;                      ## [m/s/s] acceleration of gravity
  $X::mm = .0289644;                    ## [kg/mole] molar mass of dry air
  $X::Rgas = 8.31432;                   ## [J/K/mole] gas constant
  my $inch = 0.0254;                    ## [m] definition of inch
  $X::foot = 12 * $inch;                ## [m]
  $X::gram = .001;                      ## [kg]

  $X::freezing = 273.15;                ## [K] centigrade - kelvin offset
  $X::T0 = 15 + $X::freezing;           ## [K] standard sea-level temperature
#xx  $X::T0 = 15.00520698802 + $X::freezing;    ## [K] standard sea-level temperature

## The following three quantities are believed to
## be exact, by definition:
  $X::p0 = 101325.0;                    ## [pascals] standard sea-level pressure
                                        ##  == 1013.25 mb
  $X::p0e = 760.0 / 1000 / $inch;       ## [inHg] standard sea-level pressure
                                        ##  == approx   29.92126
  $X::inHg = $X::p0 / $X::p0e;          ## [pascals]
##  NIST: quote: P/in Hg = 0.2952999 × P/kPa
#    which is consistent within roundoff error to
#    the value I calculated from 760 torr, which is considered exact
##  my $NIST_inHg = 1000/0.2952999;     ## [Pa]
##  my $delta = $X::inHg - $NIST_inHg;
##  my $inv = 1000/$NIST_inHg;
##  print "inHg: $X::inHg,  NIST: $NIST_inHg, \n  delta: $delta,  inv: $inv\n";

## The T group is valid in the troposphere only:
  $T::lapse = .0065;                    ## [K/m] standard lapse rate (troposphere)
#xx  $T::lapse = 0.00650125595831986;
  $T::nn = $T::lapse * $X::Rgas / $X::g / $X::mm;       ## exponent in station-pressure eqn
  my $rat =  $T::nn / 0.1903;
  $T::kk = $X::foot*$T::lapse/$X::T0*$X::p0e**$T::nn;   ## scale factor in station-pressure eqn
  $rat = $T::kk / 1.313e-5;
  #print " kk: $T::kk, 1.313e-5, rat: $rat\n";
  my $scale_height = $X::T0 / $T::lapse / $X::foot;
  my $sc2 = $X::p0e**$T::nn / $T::kk ;
  my $sc_faa = $X::p0e**0.1903 / 1.313e-5;
  if ($X::verbose) {
    print "nn: $T::nn,  nn': 0.1903,  rat: $rat\n";
    print "scale: $scale_height, sc2: $sc2, sc_faa: $sc_faa  nn: $T::nn\n";
    print "p0e: $X::p0e \n";
  }

  my $lam = $T::nn * $X::g * $X::mm / $X::Rgas;
  my $lam_faa = .1903 * $X::g * $X::mm / $X::Rgas;
  #print "lam: $lam,  lam_faa: $lam_faa\n";

}

sub baro {

## OK, let's work on today's instance of the problem.
## Pick up variables from the command line:

## TO DO: If we want to get fancy, this program
## could go fetch the METAR and parse it, to
## save the user the effort of doing so.

  if ($ARGV[0] ||''  =~ m'^-v') {
    $X::verbose++;
    shift @ARGV;
  }

  if (!defined($ARGV[0]) || $ARGV[0] =~ m'^-h') {
    if (!$X::verbose) {
      print $baro_msg;
    } else {
      print $baro_msg;
      print $baro_details;
    }
    exit(0);
  }

  my $Hloc = $ARGV[0] || 0;             ## local elevation in feet
  my $aset = $ARGV[1] || 29.92126;      ## altimeter setting
  if ($aset > 300) {$aset /= 100}       ## shorthand: convert 2992 to 29.92
  my $Tst = $ARGV[2];                   ## [C] station temperature
  if (!defined($Tst)) {$Tst = 15;}
  $Tst += $X::freezing;                 ## [K]
  my $Hst = $ARGV[3];                   ## [ft] station elevation
  if (!defined($Hst)) {$Hst = $Hloc;}

## My general policy: work in SI.
## Convert everything to SI at the first opportunity:
  $Hloc *= $X::foot;                    ## elevation now in meters
  $Hst *= $X::foot;
  $aset *= $X::inHg;                    ## setting now in pascals

## Suppose we are above the station.
## When the air is colder than normal,
## the local pressure will be lower than normal,
## i.e. the pressure altitude will be higher than normal:
  my $Heff = $Hst +  ($Hloc - $Hst) * ($X::T0 - $T::lapse * $Hst) / $Tst;

## TO DO:  If we want to get fancy, we should correct
## for humidity also.  Reference:
## [1] http://wahiduddin.net/calc/density_altitude.htm
## apparently based on
## [1b] http://www.meted.ucar.edu/export/asos
##
## Of course we always like:
## [2] http://williams.best.vwh.net/avform.htm
##
## Additional references, with a slight variant of the
## aset formula:
## [3] http://www.aprweather.com/pages/calc.htm
## [4] http://www.srh.noaa.gov/elp/wxcalc/altpresssc.html

## The only differences in [3] and [4] relative to [1] are:
## (a) they use 288 instead of 288.15 as the standard sea-level
##  temperature, which introduces an error less than 0.05%,
##  which I'm not gonna worry about, and
## (b) they shift the station pressure by 30 Pa (0.3 mb) for some
##  unexplained reason.  This introduces another error less than
##  0.05% (in the opposite direction).  I'm not gonna worry about
##  this, either.

## It is not important for the ASOS to use a theoretically-
## elegant formula;  it is important for the ASOS to match the
## behavior of the existing altimeters.  In turn, our calculation
## should match whatever the ASOS is doing.

## According to reference [1], this is the formula used
## by the ASOS machine:
  my $ploc = $X::p0 * (($aset/$X::p0) ** $T::nn - $Heff * $T::lapse / $X::T0) ** (1/$T::nn);
  printf ("%d pascal", $ploc);

  my $ref = $ENV{'ref'};
  if (defined($ref)){
    my $delta = $ploc - $ref;
    my $pct = 100 * $delta / $ref;
    printf ("  ... off by %d pascal == %5.2f%%", $delta, $pct);
  }
  print "\n";


};      ######## end of main program


sub HALT {
  stack_qnh($X::T0);
  stack_qnh($X::T0 - 20);
}


# Show what the scale height would be, if the troposphere were isothermal:
sub scale {
    my $mm =  $X::mm;           # for air
#    $mm = 44.01 * $X::gram;    ## for CO2
    $mm = 16.01 * $X::gram;     ## for CH4
    my $foo = $X::Rgas * $X::T0 / $X::g / $mm;
    my $bar = $foo / $X::foot;
    printf "Isothermal scale height, base e: %d m == %d ft\n", $foo, $bar;
    $foo *= log(2);
    $bar *= log(2);
    printf "Isothermal scale height, base 2: %d m == %d ft\n", $foo, $bar;
    exit(0);
}

sub stack_qnh {
  my ($tref) = @_;
  print "tref: $tref\n";
  my $Psl = P_layer(0, 0, $X::p0, $tref, $T::lapse);
  my $fmt = "%5.3f";
  $fmt = "%5.2f";
  $Psl = sprintf($fmt, $Psl / $X::inHg);
  for my $hgt (0, 2500, 5000, 7500, 10000) {
    my $press = P_layer($hgt*$X::foot, 0, $X::p0, $tref, $T::lapse);
    $press /= $X::inHg;
    my $qnh = qnh($hgt, $press);
    my $qnha = qnh_awos($hgt, $press);
    my $qnh2 = sprintf("%5.2f", $qnh);
##    $qnh2 = $qnh;
    my $pressalt = alt_indic($press, $X::p0e) / $X::foot;
    my $hind = alt_indic($press, $qnh2) / $X::foot;
    my $hxx = alt_indic($press, $qnh) / $X::foot;
    my $hxxa = alt_indic($press, $qnha) / $X::foot;
    $qnh = sprintf($fmt, $qnh);
    my $delta = sprintf("%5.3f", $hxx-$hxxa);
    my $hgtx = sprintf("%5d", $hgt);
    $press = sprintf("%5.2f", $press);
    $pressalt = sprintf("%5d", $pressalt);
    $hind = sprintf("%5d", $hind);
    $delta = "";        ###### defeat for publication
    print "A: $hgtx  Aind: $hind $delta Apr: $pressalt  ",
        "P: $press  Psl: $Psl  Qnh: $qnh \n";
  }
}

## Pressure within a layer, as a function of height.
## Physics model:  standard or nonstandard atmosphere,
##  depending on what parameters you pass in.
## Heights in meters, pressures in Pa.
## As always, $lambda is positive in the troposphere,
##  and zero in the first part of the stratosphere.
##
sub P_layer {
  my ($hh, $hb, $Pb, $Tb, $lambda) = @_;
  if ($lambda) {
##  (1/N) &:=& {g\,m_M \over \lambda\,R}
    my $N = $lambda * $X::Rgas / $X::mm / $X::g;
    return $Pb * (($Tb - $lambda*($hh - $hb)) / $Tb) ** (1/$N);
  } else {
    return $Pb * exp(-$X::g * $X::mm / $X::Rgas / $Tb * ($hh - $hb));
  }
}


## 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.
sub T_layer {
  my ($hh, $hb, $Pb, $Tb, $lambda) = @_;
  return $Tb - $lambda*($hh - $hb);
}

## Make a nice little table of density and cas/tas corrections
sub do_sweep {
  my ($hmax) = @_;
  $hmax //= 65e3;
  my $rho0 = $X::mm * $X::p0 / $X::T0 / $X::Rgas;
  print '   h/ft &   T/C &     P/Pa   &    P/inHg & ρ/(kg/m3) &      ρ/ρ0 &       cas &       tas \\\\', "\n";
  for (my $hgtft = -1000; $hgtft <= $hmax; $hgtft += 1000) {
    my ($press,$temp) = PhTh($hgtft * $X::foot);
    my $rho = $X::mm * $press / $temp / $X::Rgas;

    my $hx = sprintf("%5d", $hgtft);
    my $tx = sprintf("%5.1f", $temp - $X::freezing);
    my $fmt = "%8.0f  ";
    $fmt = "  %8.1f" if $press<10000;
    $fmt = "  %8.2f" if $press<1000;
    $fmt = "  %8.3f" if $press<100;
    my $px = sprintf($fmt, $press);
    my $p_inx = sprintf("%9.2f", $press / $X::inHg);
    my $rhox = sprintf("%9.4f", $rho);
    my $rho_ratx = sprintf("%9.4f", $rho / $rho0);
    my $cas_tasx = sprintf("%9.4f", sqrt($rho / $rho0));
    my $tas_casx = sprintf("%9.4f", sqrt($rho0 / $rho));
    print  "  $hx \& $tx \& $px \& $p_inx \& $rhox \& $rho_ratx \&"
                 . " $cas_tasx \& $tas_casx \\\\\n";
  }
}

## Make a nice little table of pressure versus pressure_altitude
## by calling the PhTh() function.
sub tabulate_Ph {
  for (my $hgt = -1000; $hgt <= 32000;) {
    my ($press,$temp) = PhTh($hgt);

    if (1) {
      my $hgtx = sprintf("%5d", $hgt);
      my $fmt = $press>10000 ? "%9.2f " : " %9.3f";
      $press = sprintf($fmt, $press);
      print  " { $press, $hgtx },\n";
    }
    if (0) {
      my $hgtx = sprintf("%9.2f", $hgt / $X::foot);
      $press /= $X::inHg;
      my $diff = sprintf("%5.2f", 29.92 - $press);
      $press = sprintf("%7.4f", $press);
      ###print  "  $diff, $press, $hgtx\n";
      print  " { $press, $hgtx },\n";
    }
    if ($hgt < 6000) {
      $hgt += 500;
    } elsif ($hgt < 11000) {
      $hgt += 1000;
    } else {
      $hgt += 3000;
    }
  }
}

## similar to above, but formatted for the brain-dead table
## in FGFS Environment/environment.cxx
sub interpolate_Ph_fgfs {
  foreach my $hgt (
    -3000,
    0.00, 2952.76, 5905.51, 8858.27, 11811.02, 14763.78, 17716.54,
    20669.29, 23622.05, 26574.80, 29527.56, 32480.31, 35433.07, 38385.83,
    41338.58, 44291.34, 47244.09, 50196.85, 53149.61, 56102.36, 59055.12,
    62007.87,
    65000, 68000, 71000, 74000, 77000, 80000, 83000, 86000, 89000,
    92000, 95000, 98000, 101000) {
    my ($press,$temp) = PhTh($hgt*$X::foot);

    if (1) {
      my $hgtx  = sprintf("%9.2f", $hgt);
      my $presx = sprintf("%7.4f", $press/$X::p0);
      my $tx    = sprintf("%7.3f", $temp/$X::T0);
      print  " { $hgtx, $tx, $presx },\n";
    }
    if (0) {
      my $hgtx = sprintf("%9.2f", $hgt / $X::foot);
      $press /= $X::inHg;
      my $diff = sprintf("%5.2f", 29.92 - $press);
      $press = sprintf("%7.4f", $press);
      ###print  "  $diff, $press, $hgtx\n";
      print  " { $press, $hgtx },\n";
    }
  }
}


## Altimeter setting.
## height in feet
## pressure in inHg
## Troposphere only
sub qnh{
  my ($h, $P) = @_;

#   AS = (Pa^N + K Ha)^(1/N)
    my $nn = $T::nn;
    my $kk = $T::kk;

##  return $P * (1 + $kk * $h / $P **  $nn) ** (1/$nn);
    return ($P ** $nn + $kk * $h) ** (1/$nn);
}

sub qnh_awos{
  my ($h, $P) = @_;

#   AS = (Pa^N + K Ha)^(1/N)
    my $nn = $T::nn;
    my $kk = $T::kk;
    $nn = .1903;
    $kk = 1.313e-5;

##  return $P * (1 + $kk * $h / $P **  $nn) ** (1/$nn);
    return ($P ** $nn + $kk * $h) ** (1/$nn);
}

sub leftovers {

  my $h = 5000;
  print "Pressure at $h: ", P_layer($h, 0,     $X::p0, $X::T0, $T::lapse), "\n";

  $h = 20000;
  print "Pressure at $h: ", P_layer($h, 11000, 22632.1, 216.65, 0), "\n";

  $h = 86000;
  my $ppp = P_layer($h, 71000, 3.95642, 214.65, 0.002);
  my $foo = $ppp/$X::inHg;
  print "Pressure at $h:  $ppp / $foo \n";
  print "Temperature: ", 214.65 -0.002*(86000-71000) . "\n";
  my $qnh = qnh(-1000, 29.92);
  $qnh = sprintf("%1.2f", $qnh);
  printf ".... qnh: $qnh\n";

}

## Indicated Altitude
## pressure in inHg
## qnh in inHg
## Troposphere only
sub alt_indic {
  my ($p, $qnh) = @_;
  return $X::T0 * (($qnh/$X::p0e)**$T::nn - ($p/$X::p0e)**$T::nn) / $T::lapse;
}

##################################################################333
# leftovers

## Pressure as a function of height.
## Valid below 32000 meters,
##   i.e. troposphere and first two layers of stratosphere.
## $hh in meters, pressure in Pa.
sub Ph_older {
  my ($hh) = @_;
  if ($hh <= 11000) {
    return P_layer($hh, 0, $X::p0, $X::T0, $T::lapse);
  } elsif ($hh <= 20000) {
    return P_layer($hh, 11000, 22632.06, 216.65, 0);
  } elsif ($hh <= 32000) {
    return P_layer($hh, 20000,  5474.89, 216.65, -0.001);
  }
  return -9e99;
}

## Pressure as a function of height.
## heights in meters, pressures in Pa.
sub Ph_old {
  my ($hh) = @_;
  my $maxrow = $#X::ISA_table;
  for (my $ii = 0; $ii <= $maxrow; $ii++){
    my $row = $X::ISA_table[$ii];
    my ($ndx, $hgt, $hgt_ft, $p0, $p0_in, $t0, $lapse, $lapse_ft) = @$row;
    my $xhgt;
    if ($ii < $maxrow) {
      my $xrow = $X::ISA_table[1+$ii];
      $xhgt = $xrow->[1];
    }
    if (!defined $xhgt || $hh <= $xhgt ) {
      return P_layer($hh, $hgt, $p0, $t0, $lapse);
    }
  }
  die "Ran out of layers";    # should never happen
}

## Temperature as a function of height.
## heights in meters, pressures in Pa.
sub Th_old {
  my ($hh) = @_;
  my $maxrow = $#X::ISA_table;
  for (my $ii = 0; $ii <= $maxrow; $ii++){
    my $row = $X::ISA_table[$ii];
    my ($ndx, $hgt, $hgt_ft, $p0, $p0_in, $t0, $lapse, $lapse_ft) = @$row;
    my $xhgt;
    if ($ii < $maxrow) {
      my $xrow = $X::ISA_table[1+$ii];
      $xhgt = $xrow->[1];
    }
    if (!defined $xhgt || $hh <= $xhgt ) {
      return T_layer($hh, $hgt, $p0, $t0, $lapse);
    }
  }
  die "Ran out of layers";    # should never happen
}

@X::ISA_table_wiki_orig = (
#vvv  ............... boundary ......................... vvvv layer above bdy vvvvvvv
# 0    1      2          3          4             5        6            7
#(b)  (m)    (ft)       (pascals)  (inHg)        (K)      (K/m)        (K/ft)
[0,   0,       0,         101325,  29.92126,     288.15,  -0.0065,     -0.0019812   ],
[1,   11000,  36089,      22632,   6.683245,     216.65,  0.0,         0.0          ],
[2,   20000,  65617,      5474,    1.616734,     216.65,  0.001,       0.0003048    ],
[3,   32000,  104987,     868,     0.2563258,    228.65,  0.0028,      0.00085344   ],
[4,   47000,  154199,     110,     0.0327506,    270.65,  0.0,         0.0          ],
[5,   51000,  167323,     66,      0.01976704,   270.65,  -0.0028,     -0.00085344  ],
[6,   71000,  232940,     4,       0.00116833,   214.65,  -0.002,      -0.0006097   ]
);