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

use strict;

my $usage_msg =  <<EoF;
Program to calculate various atmosphere-related and
barometer-related information.

Usage:
  barometry -option ...

Option examples:
  -help              # print this message
  -baro              # find local barometric pressure
  -baro -h           # info about baro option
  -baro -v -h        # details about baro option
  -formula           # explain the main formula used
  -scale             # scale height for an isothermal atmo
  -halt              # example of high altimeter @ low temp
  -cas               # tabulate tas and cas versus altitude
  -tas               # tabulate tas and cas versus altitude
  -check1 p0 t0      # check the PhTh function
  -check2            # print the table that defines ISA
  -tex               # ISA table in latex format
  -cxx               # ISA table in c++ format
  -tab               # h as a function of p, c++ format
  -interp            # p and t as a function of h, c++
EoF

my $baro_msg = <<EoF;
Calculate the local barometric pressure [in pascals]
based on nearby weather reports.  This is useful if you
want to know the pressure but don't own a precision
barometer, and don't want to buy one.  Even if you
do own a barometer, it's nice to have an independent
cross-check.

Usage:
  barometry -baro Hloc aset Tst Hst
Arguments are:
     Hloc       your local elevation [feet]
     aset       the reported altimeter setting [inches Hg]
     Tst        the reported temperature at the station [C]
     Hst        the elevation of the station [feet]

The Tst and Hst arguments are optional.  They improve
accuracy, especially in very hot or very cold weather,
especially if your local elevation differs from the
reporting station elevation.

If you don't know your local elevation, you can figure
it out with the help of a topo map.  Don't forget to
account for the height of your workspace relative to
ground level.

You can also get rough elevation info from GPS, but
beware that its vertical uncertainty is at least ten
times worse than its horizontal uncertainty.  Averaging
for a long time might help.

You can get the altimeter setting from a METAR report.
See below for data sources.

Usage example:
   # Pressure at the Prescott airport,
   # computed from observations made right there:
   # METAR KPRC 201453Z 00000KT 10SM CLR 23/M02 A3020
   barometry 5045 30.20
   84977 pascal

   # Same thing,
   # inferred from observations made at Flagstaff
   # METAR KFLG 201456Z 06005KT 10SM CLR 22/00 A3037
   ref=84977 barometry 5045 30.37 22 7014
   85029 pascal  ... off by 52 pascal ==  0.06%

   # And the other way around: the pressure at Flagstaff
   # computed from the local observation there:
   barometry 7014 30.37
   79375 pascal

   # Same thing, inferred from Prescott observations:
   ref=79375 barometry 7014 30.20 23 5045
   79267 pascal  ... off by -107 pascal == -0.14%

Note that this is a pretty strict test, since the two
locations are 50 miles apart, and differ by 2000 feet of
elevation.  Most of the population should be able to get
a METAR from somewhere not nearly so far away.

In normal usage, as opposed to a test, you would just
calculate the pressure at your location (not multiple
locations) by feeding this program the altimeter
settings from some nearby reporting stations, and
averaging the results.

In addition to this perl version, there is also a
spreadsheet version:
  http://www.av8n.com/physics/barometry.xls

=====================
Data sources:

1a) Super-convenient, bookmarkable URL for your favorite
observation site, worldwide (not just US):
  http://weather.noaa.gov/weather/current/KFLG.html
  http://weather.noaa.gov/weather/current/KPRC.html
  http://weather.noaa.gov/weather/current/EGLL.html

(just plug in four-letter ICAO identifier for the site you
want, in the URL at the place indicated).

1b) One-step access to worldwide METAR data (just the
METAR, without much other stuff)
  http://weather.noaa.gov/weather/metar.shtml
Example:  Enter EGLL to get the METAR for Heathrow.

1c) Two- or three-step access to US data (less convenient
if you already know the ICAO code, but more convenient
otherwise):
  http://weather.noaa.gov/index.html
Just enter the name of your state in the first box and
go from there.

This will lead you to a nice bookmarkable URL of the type
mentioned in item (1a).

2) One-step access to history (going back two days) of METAR
observations for a couple thousand sites ... mostly US,
Canada, and Caribbean:
  http://www.srh.noaa.gov/data/obhistory/

(Results are bookmarkable in the obvious way).

3) Key for decoding METARs:
  http://www.met.tamu.edu/class/METAR/quick-metar.html

4) A comprehensive list of weather reporting sites
(including METAR sites of course) is at
  http://adds.aviationweather.gov/metars/stations.txt
This includes the conventional place-name, ICAO
identifier, lat/long, and elevation (in meters).

5) A map of weather stations with phone numbers:
  http://www.faa.gov/asos/map/map.htm
You can call to get an observation that is current
up-to-the-minute (as opposed to online data, which
is generally updated only once per hour, a few minutes
before the top of the hour).

EoF

my $formula = <<EoF;
I asked the folks who manufacture AWOS machines to tell
me what algorithm they use for computing Altimeter Setting.

They replied:
> 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 = <<EoF;
===========
Details:

I figure the vast majority of the US population lives within 30 miles
of a METAR reporting station ... probably more than one.

My advice: Rather than theorizing too much about the details of the
location-dependence of the model, just use this program, as described
above.  Apply it to multiple nearby stations and check for
consistency.  If the results are not too scattered, and the stations
are not wildly far away (horizontally or vertically), you should
should be able to take some kind of average without worrying about the
details.


However, for those who dote on details:

2) The second part of the task is to convert field pressure
to your local pressure.  To do this, my code uses a model
that may be described, roughly, as being second-order in
altitude and zeroth-order in latitude and longitude.

  BTW, when I refer to "my code", I mean both this perl
  script and the corresponding spreadsheet:
    http://www.av8n.com/physics/barometry.xls

2a) Considering only altitude for a moment: The zeroth-order
model would estimate your local pressure to be the same as
the station pressure.  The first-order term corrects for the
difference in elevation between the station and your locale,
basing this correction on the standard atmosphere.  As a
correction to this correction, we account for non-standard
temperature.  At this point we have a second-order
description, in the sense that it depends on the product of
two smallish quantities: the difference in elevation, and
the discrepancy between temperature and standard
temperature.

If you want to take into account nonlinearities in the
temperature profile (inversions and such) you can do that,
I suppose, but that's getting into third-order territory.
My code doesn't do that.


2b) Let's now consider lattitude and longitude. The
zeroth-order model is that your local pressure should be
independent of the lat/long of the station(s) from which you
obtained the raw data.

If you use multiple stations (as you should, if you
can), the spreadsheet averages them ... which means it
is applying the zeroth-order model.  The perl script
just processes a single station;  you get to do the
averaging yourself.

Given enough stations and information about their
location, you could implement some sort of interpolation
scheme.  Possible methods in the ad-hoc category are:
 -- perhaps some sort of polynomial fit, or
 -- perhaps some sort of Parzen windows.
These are way beyond what my code does.

Even better would be to use a non-ad-hoc method, based on
real physics and real data.  In the atmosphere, pressure
gradients are related to wind velocity.  There's lots of
data available about regional-scale wind and pressure
patterns.  So you could do all sorts of fancy modeling.

On the other hand, if your laboratory sits atop an
airfoil-shaped hill, there will be pressure variations that
depend on wind.  These will be sensitive to the exact shape
of the hill, in ways that aren't reflected in the regional
wind/pressure patterns.  Indeed, if there are gusty winds,
there will be pressure fluctuations that invalidate the
assumption that your experiments are carried out under
constant-pressure conditions.  For wind speeds less than 45
mph, these effects are unlikely to exceed 250 pascals (and
for 20 mph, 50 pascals ... it's quadratic).  You can decide
for yourself whether this is significant for your
application.

EoF


main: {                 ## set scope of local-to-main variables

  $X::verbose = 0;

  argloop: for (;;) {

    my $arg = shift @ARGV || '-help';

    if ($arg eq '-hh' || $arg eq '-H') {
      print $baro_msg;
      print $baro_details;
      exit(0)
    }

    if ($arg =~ m'^-h') {
      print $usage_msg;
      exit(0);
    }

    if ($arg eq '-v') {
      $X::verbose++;
      next argloop;
    }

    set_nature();

    if ($arg eq '-formula') {
      print $formula;
      exit(0);
    }

    if ($arg eq '-scale') {
      scale();
      exit(0);
    }


    if ($arg eq '-halt' || $arg eq '-HALT') {
      HALT();
      exit(0);
    }

    if ($arg eq '-check' || $arg eq '-check1') {
      my $arg0 = $ARGV[0] ? $ARGV[0] * $X::inHg : undef();
      my $arg1 = defined $ARGV[1] ? $ARGV[1] + $X::freezing : undef();
      check_PhTh($arg0, $arg1);
      exit(0);
    }

    if ($arg eq '-tab' ) {
      tabulate_Ph();
      exit(0);
    }

    if ($arg eq '-tex' ) {
      tex();
      exit(0);
    }

    if ($arg eq '-cxx' ) {
      cxx();
      exit(0);
    }

    if ($arg eq '-check2' ) {
      check2();
      exit(0);
    }

    if ($arg =~ m'^-cas' || $arg =~ m'^-tas') {
      tabulate_rho_cas_tas();
      exit(0);
    }

    if ($arg =~ m'^-baro') {
      baro();
      exit(0);
    }

    if ($arg =~ m'^-interp') {
      interpolate_Ph_fgfs();         ## for FlightGear
      exit(0);
    }

    die "Extraneous verbiage: '$arg'\n";

  }

  exit(0);

}

## Pressure as a function of height.
## 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 tabulate_rho_cas_tas {
  my $rho0 = $X::mm * $X::p0 / $X::T0 / $X::Rgas;
  for (my $hgtft = -1000; $hgtft <= 36000; $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 $px = sprintf($press>10000 ? "%9.0f " : " %9.0f", $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, 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);
      $press = sprintf($press>10000 ? "%9.2f " : " %9.3f", $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   ]
);