#! /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: 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; 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, apply this program to data from 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_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_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 ] );