| MTP Home Page |

How to Determine Whether In-Situ Temperature and Pressure Measurements are Accurate

MJ Mahoney

Last Revision: March 5, 2008

Abstract

It is often assumed that in-situ temperature and pressure measurements aboard atmospheric research aircraft are accurate when in fact they are not. This is often the case when PI-led temperature and pressure measurements are not available, because no one is responsible for maintianing the measurement systems. Nearly a decade ago we introduced a simple technique to validate the temperature and pressure measurements made by the Navigation Data Recorder (NDR) aboard the NASA ER-2. It involves comparing the difference between geometric (Zg) and pressure (Zp) altitudes (as a function of either geometric or pressure altitude) based on the in-situ aircraft measurements to the same quantity (Zg-Zp) for a radiosonde launched near the aircraft's flight track.

Outline

Introduction

There are three flavors of altitude that will be involved in this discussion: pressure altitude, geopotential height, and geometric altitude. The pressure altitude scale is based on the International Civil Aviation Organization's (ICAO) International Standard Atmosphere (ISA), which is the same as the US Standard Atmosphere (1976) to an altitude of 32 km. It is based on an average mid-latitude temperature profile represented by piece-wise continuous linear segments between different geopotential heights. Geopotential heights can also be calculated from radiosonde pressure, temperature and humidity measurements using the hypsometric equation (see Ref. 1, Equation 7). The difference between pressure altitude and geopotential height is that the former is defined by an average mid-latitude standard atmosphere (that is, not the real one that an aircraft is flying in, but nevertheless an agreed upon standard), and that standard atmosphere has zero humidity.

How are pressure altitude and geopotential height related to geometic altitude? This actually gets a little complicated, and is discussed in detail elsewhere. Historically, mean sea level (MSL) has been used as the zero of geometric altitude (or elevation), and this is closely approximated by an equipotential surface for the earth's gravity field called the geoid. Although much smoother than the topographic surface of the earth, the geoid has significant vertical undulations due to the large scale distribution of Earth's mass. A much smoother reference surface is the reference ellipsoid, which is used to the approximate the shape of the earth. The reference ellipsoid is convenient because three dimensional coordinates (latitude, longitude, and altitude) are easily defined with respect to it. For example, the Global Positioning System (GPS) of satellites, which measures geometric altitude, uses the WGS-84 ( World Geodetic System 1984 ) reference ellipsoid. Time corrections are transmitted via the GPS satellites, which work with the WGS-84 model software inside the GPS receivers. GPS measurements made relative to this reference ellipsoid are generally converted to different vertical datums by the software internal to the GPS receiver. This involves a geoid model, the current one being EGM96 ( Earth Gravitational Model 1996 ). The GPS receiver allows users to select between many regional vertical datums as well; that is, different vertical reference surfaces.

To evaluate the accuracy of in situ aircraft temperature and pressure measurements, we need to be able to calculate Zg-Zp using a nearby radiosonde's measurements and also the aircraft measurements. This is trivial for the aircraft measurements because both the geometric altitude (GPS altitude) and pressure altitude are readily available. It is a little more involved for the radiosonde measurements. Pressure altitude is readily calculated from the radionsonde pressure measurements, and geopotential height is readily calculated from the pressure, temperature and humidity measurements. What remains to be done is to convert geopotential height to geometric altitude.

Converting Between Geopotential Height and Geometric Altitude

As discussed in more detail elsewhere, this can be done by writing the expression for the geopotential height between two levels (H12):
(1) GeopotentialHeight.png
in differential form:

(2) Gamma Gk 45 dh = Gamma Gk (z,phi) dz

In fact, this expression is often used as the definition of the geopotential height. Denoting for simplicity R as the radius of the earth at latitudephi, and using the inverse square law for gravity, we can write for any latitude (phi ) and geopotential height (z):

(3) Gamma Gk (R + z,phi) = Gamma Gk (R,phi ) (R/ (R + z))2.

Substituting equation (3) into equation (2) and integrating from z = 0 to z = Z -- assuming Gamma Gk (z,phi ) does not vary with z -- gives:

(4) H_Z_PHI_R

Now of course Gamma Gk (z,phi ) does vary with geometric altitude (z) so this must be dealt with. A trick suggested by W. D. Lambert in 1949 in the Smithsonian Meteorological Tables (SMT), which is a  frequently cited source for converting between geometric altitude and geopotential height, was to compensate for the variation of Gamma Gk (z,phi ) by letting R assume an appropriate value that satisfied certain boundary conditions. This would account for the fact that equation (3) only applies to a non-rotating sphere composed of spherical shells of equal density, which is not true for Earth, and for the fact that the centrifugal component of gravity increases linearly with radius, rather than inversely with the square of the radius as shown in this equation. Taking the derivative of equation (3) with respect to z, evaluating the resulting expression at z = 0, and then solving for R, we obtain:

(5) R_SAO

Note the minus sign in the denominator of equation (5), and remember that the derivative is evaluated at z=0 (which is why the z dependence is dropped on Gamma GkSMT). The notation RSMT is used because this technique was originally described in the Smithsonian Meteorological Tables. Note that RSMT is not the radius of the Earth at any particular latitude (as is often mistakenly assumed to be the case); rather it is a value that is needed to account for the combined effects of gravitational and centrifugal forces with geometric altitude z. Using this notation, equation (4) becomes:

(6) H_SMT_Z_PHI ,

or on inverting to solve for Z

(7) Z_SMT_H_PHI .

Visual Basic code to perform these calculations is provided below, along with code to calculate pressure altitude from pressure or vice versa. So the reader has everything that is needed to proceed.

Behaviour of Zg-Zp

As shown in Figure 1a and 1b the difference between the geometric altitude (Zg) and pressure altitude (Zp) will take different shapes as a function of either Zg or Zp depending upon the temperature structure in the atmosphere. Three temperature profiles are shown in Figure 1a: a tropical one, a mid-latitude one (for which we use the US Standard Atmosphere 1976), and an Arctic winter one. 

ModelTprofile.png
Figure 1a. Representative tropical (blue), mid-latitude (red) and polar (green) temperature profiles.
 Zg_Zp_Model.png
Figure 1b. The difference between geometric and pressure altitude as a function of geometric altitude for the three representative temperature profiles in Figure 1a.

We note in Figure 1b that the difference between geometric and pressure altitude (Zg-Zp) is smallest for the mid-latitude temperature profile, which is represented by the US Standard Atmosphere (1976).

Consider for the moment, the Arctic winter profile (green). By way of example, when the NASA ER-2 is flying at a pressure altitude of 20 km (~54 hPa) in the Arctic winter, this pressure level will be at a much lower geometric altitude than if the ER-2 were flying at the same pressure level in the tropics simply because cold Arctic atmosphere is much more compressed that the warm tropical atmosphere. So Zg-Zp is negative in a cold atmosphere and positive in a warm atmosphere. At mid-latitudes Zg-Zp is smallest because the geopotential height was defined to agree closely with pressure altitude. Only the small latitude and altitude dependencies of gravity come into effect when calculating Zg.

The temperature and pressure altitude reported by the ER-2 Nav Data Recorder (NDR) during TC4 is unuseable for science. Indeed, issues with the NDR were first noticed during the first SOLVE capaign in 2000, and then again during the CAMEX-4 campaign in 2001. The problems have become worse with time. We were able to apply constant pressure corrections during the earlier campaigns but this technique no longer works.

Zg-Zp_20010909.GIF
Figure 2. A plot of  Zg-Zp versus Zp for the CAMEX-4 flight of September 9, 2001.

Figure 2 shows an example from the CAMEX-4 campaign of how the Zg-Zp verification technique works. (Figure 2 which was analyzed seven years ago, plots Zg-Zp vs Zp. Today we more commonly plot Zg-Zp vs Zg, so that the potentially problematic measurement, pressure altitude, appears on only one axis.)  Shown in dark and light blue are the Zg-Zp for two radiosondes launched from Key West, FL, before and after the ER-2 flight. For these radiosondes, Zg is calculated by applying a latitude-dependent correction to the reported geopotential height using Equation (7). These radiosondes are the standard to which the data from the ER-2 Nav Data Recorder (NDR) are to be compared. There is also a green trace showing for the DC-8 derived Zg-Zp. The dark brown trace is the difference between the ER-2 GPS and pressure altitude (Zg-Zp) as a function of the pressure altitude (Zp), and it does not agree with the radiosonde and DC-8 data. By applying a small correction (-2.3 hPa) to the reported ER-2 pressure altitude, we find that the corrected ER-2 Zg-Zp agrees quite well (orange trace) with the two radiosondes above 9 km. The jump at 9 km is caused when the ER-2 Environmental Control System switches from isobaric to differential control mode; this is discussed in the next section. The exact altitude that this occurs at can vary depending on the sequence in which the pressure is adjusted in each payload area.

A Model for the NDR Zp Error
ECS
Figure 3. ER-2 Environmental Control Sytem (ECS) Pressure Altitude Profile

As is shown in Figure 3, under normal operation the ER-2 Environmental Control System (ECS) will allow the cockpit and payload areas to remain unpressurized below 7900 feet (2.4 km). Above this altitude the ECS switches to an isobaric control mode in which a pressure altitude of 7500 feet (2.3 km) is maintained until 18,300 feet (5.6 km). Above 5.6 km, a differential control mode is maintained in which the cabin/payload areas is kept 3.88 psi above the ambient outside air pressure. This is where things get a little unpredicable, because there is separate valving for the cockpit and each payload area, and (needlesstosay) the cockpit always takes priority. Under differential control mode , the cockpit will be at a pressure altitude of 28,381 feet (8.647 km, 4.7 psi, 323.9 hPa) when the outside air pressure altitude is 65,000 feet (19.805 km, 0.82 psi, 56.46 hPa); that is, a 3.88 psi difference.
PiPoPp
The temporal and Zp dependence of the pressure altitude correction can (at least qualitatively) be understood by the following simple model for a leak in the static pressure line (see the figure to the right). Since we are interested in understanding the error at cruise altitudes, we can assume that the ECS is in differential pressure control mode . If the ECS is working perfectly (that is, no delays), then:

Pi - Po = 3.88 psi = 267.44 mb

Therefore, the effect of a leak should be constant with time, since the pressure gradient through the leak is constant! So in this perfect scenario, there should be a constant positive pressure offset on at the pressure transducer (Pp) at all altitudes above 5.6 km, the altitude where the handoff to differential control mode takes place. When best fits are made to all the CAMEX-4 flights, the increased pressure that the pressure transducer sees is: Pcorr = 2.28 +/- 0.17 hPa (see this web page for the CAMEX-4 analysis); that is, Pp = Po + Pcorr due to the leak. It is reassuring that the variation from flight to flight is relatively small. It turns out, however, that a better fit to Zg-Zp in Figure 2 can be obtained if a term depending on the time since launch is also included. Without this term, the ascent and descent data separate from one another as discussed on the CAMEX-4 Zg-Zp analysis web page.

On ascent it will take some time for the ECS to establish the correct differential pressure, so it might be expected that the cabin pressure Pi is more than 3.88 psi greater than the outside static air pressure Po. As a result, if there is a leak in the static pressure line during ascent, the leak of cabin air Pi into the static pressure line will increase Po so that the measured pressure Pp is greater than Po + Pcorr. This should introduce a temporal term to the correction, which appears to be needed in the real data.

Conversely, on descent the cabin pressure Pi may be <3.88 psi greater outside static air pressure Po because the ECS is playing catch up. The effect of this is to reduce the effective Pcorr -- the opposite to ascent. Thus, on descent the measured Zp will be greater than the true pressure altitude. This is exactly what we see in Figure 2. Also, this delta- correction to Pcorr is generally more pronounced on descent than on ascent because descent is more rapid and the ECS is less able to keep up.

As a reminder, all of the pressure changes discussed here regarding the ECS have nothing to do with the pressure that the NDR measures. The ECS uses the pressure sensors that the pilot sees, and the Lockheed-Martin ground crew maintains these sensors. In fact, these pressure sensors were leak checked before CAMEX-4, but they are unforturnately not recorded.  The NDR sensors are not maintained, although I was recently told that following my queries after CAMEX-4 ended, these sensors were checked. There was a pressure leak found in the static air line, one of the two temperature sensors was not operational, and the calibration was invalid on the temperature probe that was working.

A further complication on the correction of the NDR pressure altitude occurs during flights when the ER-2 climbs or dips during flights, as this reverses the sense of the pressure difference between the static pressure line and the cockpit/payload areas. We have not attempted to model this effect.

The Current Situation

TC4_20070724a.png
Figure 4a. Raw ER-2 NDR data with no pressure corrections.
TC4_20070724b.png
Figure 4b. Raw ER-2 NDR data with a +1.2 hPa pressure correction fixes cruise altitude data.
TC4_20070724c.png
Figure 4c. Same as Figure 4b, but with a +12 and -24 hPa  altitude dependent corrections on ascent and descent, respectively.

Unfortunately the ER-2 pressure altitude has deteriorated significantly since the CAMEX-4 campaign. It is no longer possible to apply a constant pressure correction and obtain an improvement like the orange trace in Figure 2. Figures 4 shows the current situation based on data from the July 24, 2007, TC4 campaign flight. The ER-2 Zg-Zp is shown as the cyan trace in Figures 4 and the RAOB by the yellow trace. A +2.0 hPa pressure correction brought the brown ER-2 data at 20 km into agreement with the yellow MROC RAOB (Figure 4a), but the ascent data (lower cyan trace) and descent data (upper cyan trace) showed no improvement like in Figure 2. The green trace in Figures 4 provides another diagnostic.

We have also used the ER-2 pressure and temperature measurements to integrate the hypsometric equation to obtain geopotential height (just like for radiosondes) and then convert it to geometric altitude. We could not integrate this data from the known ground geopotential height because the ER-2 does not report temperatures below the freezing altitude (~5-6 km). As a result, the starting altitude adjusted to show agreement with the RAOB Zg-Zp trace at 20 km. This green trace actually agrees quite well with the yellow RAOB trace to the peak of the Zg-Zp trace. Because the erratic behaviour seen in the cyan trace does not appear in the green trace, it suggests that the erratic behaviour in fact is in the ER-2 GPS altitude measurement. This is used to calculate Zg-Zp for the cyan trace, but is not used to calculate Zg-Zp for the green trace.

The Bottom Line

The ER-2 Nav Data Recorder measurements of pressure and geometric altitude do not agree with the same data derived from radiosondes launched near the ER-2 flight track. This problem was first noticed during the first SOLVE campaign in 2000; it was still a problem during the CAMEX-4 campaign in 2001, and has degraded significantly since then. The ER-2 NDR data from the 20007 TC4 campaign is unuseable for scientific calculations.

References

1. MJ Mahoney, A Discussion of Various Measures of Altitude, 2001.
2. R. J. List, Editor, Smithsonian Meteorological Tables, 1968.

Code

Listed below is Visual Basic code that can be used calculate geometric altitude from geopotential height (and vice versa), and to calculate pressure altitude from pressure (and vice versa). It should be easily translated to other programming languages.

Convert Geopotential Height to Geometric Altitude

Function fZhtoZgSMT(Zgh!, Latitude!)
Dim r#
' Convert geopotential height (Zgh) to geometric altitude (Zg)
' Zgh ......... Geopotential height (km)
' Latitude .... Degrees
' Return ...... Geometric altitude (km)
'
' Verified against Smithsonian Meteorological Tables (page 222-223)
' if 9.80665 set to 9.8 as in tables
  r = fRsmt(Latitude)                'Ad hoc radius
  fZhtoZgSMT = r * Zgh / (fgnSMT(Latitude) * r / 9.80665 - Zgh)
End Function

Convert Geometric Altitude to Geopotential Height

Function fZgToZghSMT(Zg!, Latitude!)
' Convert geometric altitude to geopotential height
' Zg ........ Geometric altitude (km)
' Latitude .. Degrees
' Return .... Geopotential height
'
' Verified against Smithsonian Meteorological Tables (page 220-221)
' if 9.80665 set to 9.8 as in tables
Dim r!
  r = fRsmt(Latitude)
  fZgToZghSMT = r * (fgnSMT(Latitude) / 9.80665) * (Zg / (Zg + r))
End Function

Ad Hoc Radius

Function fRsmt(ByVal Latitude!) As Single
' Smithsonian Meteorological Tables p.218 Equation 6
' Ginned up earth radius to compensate for centrifugal force
' variation with latitide
' Note that this is not the Earth ellipsoid radius!
  fRsmt = -2# * fgnSMT(Latitude) / fdgdzSMT(Latitude) / 1000
End Function

Normal Gravity

Function fgnSMT(Latitude!)
' Normal Gravity vs Latitude from Smithsonian Meteorological Tables (page 488)
' rpd = radians per degree = 1.74532925199433E-02
Dim Cos2L!
  Cos2L = Cos(2 * Latitude * rpd)
  fgnSMT = 9.80616 * (1 - 0.0026373 * Cos2L + 0.0000059 * Cos2L ^ 2)
End Function

Vertical Variation in Normal Gravity

Function fdgdzSMT(Latitude!)
Dim Cos2L!, Cos4L!
' Rate of change of gravity with altitude
' From the Smithsonian Meteorological Tables p.218 Equation 7
' rpd = radians per degree = 1.74532925199433E-02
  Cos2L = Cos(2 * rpd * Latitude)
  Cos4L = Cos(4 * rpd * Latitude)
  fdgdzSMT = -(3.085462 * 10 ^ (-6) + 2.27 * 10 ^ (-9) * Cos2L - 2 * 10 ^ (-12) * Cos4L)
End Function

Convert Presssure Altitude to Pressure

Function fZtoP!(z!)
' Convert US Standard Atmosphere 1976 from Pressure Altitude [km] to Pressure [mb]
' Uses Function fTstd (z) to calculate US standard temperatures
' This approach is taken as it reduces number of calculations
' MJ Mahoney JPL 19980507
'  Z         T          P           D
' km         K         mb           kg/m3
'  0      288.15  1013.25          1.22499919057116       1
' 11      216.65   226.3206        0.363917777824827      2
' 20      216.65    54.74888       8.80347996750117E-02   3
' 32      228.65     8.680185      1.32249977610308E-02   4
' 47      270.65     1.109063      1.42753221375531E-03   5
' 51      270.65     0.6693885     8.61604682553416E-04   6
' 71      214.65     3.956419E-02  6.42109635061132E-05   7
' 84.852  186.946    3.733836E-03  6.95787870826203E-06   8
' 90      186.946 9  0.2908
' 95      188.4  10  1.34
' 100     195.1  11  2.74
' 105     208.8  12  6.24
' 110     240.0  13 12.0
' 120     360.0  14

Dim P!, T!

T = fTstd(z)
Select Case z
Case Is <= 11
' P = P(1)*(T/T(1))^(-1000*cg/(cRs*LR(1)))
' P = 1013.25 * (T / 288.15) ^ (-1000 * cg / (cRs * -6.5)
  P = 1013.25 * (T / 288.15) ^ (5.2558761507598)
Case Is <= 20
' P = 226.3206 * Exp(-1000 * cg * (z - 11.) / (cRs * 216.65))
  P = 226.3206 * Exp(-0.157688414400825 * (z - 11#))
Case Is <= 32
' P = 54.74888 * (T / 216.65) ^ (-1000 * cg / (cRs * 1.)
  P = 54.74888 * (T / 216.65) ^ (-34.1631949799387)
Case Is <= 47
' P = 8.680185 * (T / 228.65) ^ (-1000 * cg / (cRs * 2.8)
  P = 8.680185 * (T / 228.65) ^ (-12.2011410642638)
Case Is <= 51
' P = 1.109063 * Exp(-1000 * cg * (Zs(i) - Zs(i - 1)) / (cRs * 270.65))
  P = 1.109063 * Exp(-0.126226473230884 * (z - 47))
Case Is <= 71
' P = 0.6693885 * (T / 270.65) ^ (-1000 * cg / (cRs * -2.8)
  P = 0.6693885 * (T / 270.65) ^ (12.2011410642638)
Case Is <= 84.852
' P = 0.03956419 * (T / 214.65) ^ (-1000 * cg / (cRs * -2.0)
  P = 0.039564189818084 * (T / 214.65) ^ (17.0815974899694)
Case Is <= 90
' P = 0.003733834 * Exp(-1000 * cg * (z - 84.852) / (cRs * 186.946))
  P = 0.003733834 * Exp(-0.182743653140151 * (z - 84.852))
Case Is <= 95
' P = 0.001457425 * (T / 184.946) ^ (-1000 * cg / (cRs * 0.2908)
  P = 1.45742511874549E-03 * (T / 186.946) ^ (-117.480037757699)
Case Is <= 100
' P = P(11)*(T/T(11))^(-1000*cg/(cRs*LR(11)))
' P = 0.0005808211 * (T / 188.4) ^ (-1000 * cg / (cRs * 1.34)
  P = 5.8654139565495E-04 * (T / 188.4) ^ (-25.4949216268199)
Case Is <= 105
' P = 0.0001815319 * (T / 195.1) ^ (-1000 * cg / (cRs * 2.74)
  P = 2.40645796828482E-04 * (T / 195.1) ^ (-12.4683193357411)
Case Is <= 110
' P = 0.00007788813 * (T / 208.8) ^ (-1000 * cg / (cRs * 6.24)
  P = 1.03251578598705E-04 * (T / 208.8) ^ (-5.4748709903748)
Case Else
' P = 0.00003633683 * (T / 240.) ^ (-1000 * cg / (cRs * 12)
  P = 4.81695302325482E-05 * (T / 240#) ^ (-2.84693291499489)
End Select
fZtoP = P
End Function

Convert Pressure to Pressure Altitude

Function fPtoZ(ByVal P!) As Single
Dim z!
' Convert US Standard Atmosphere 1976 from Pressure [mb] to Pressure Altitude [km]
' MJ Mahoney JPL 19980507
'  Z         T          P           D                     N  LR
' km         K         mb           kg/m3                    K/km
'  0      288.15  1013.25          1.22499919057116       1  -6.5
' 11      216.65   226.3206        0.363917777824827      2   0.0
' 20      216.65    54.74888       8.80347996750117E-02   3   1.0
' 32      228.65     8.680185      1.32249977610308E-02   4   2.8
' 47      270.65     1.109063      1.42753221375531E-03   5   0.0
' 51      270.65     0.6693885     8.61604682553416E-04   6  -2.8
' 71      214.65     3.956419E-02  6.42109635061132E-05   7  -2.0
' 84.852  186.946    3.733836E-03  6.95787870826203E-06   8   0.0
' 90      186.946    1.45742511874549E-03                 9  0.2908
' 95      188.4      5.8654139565495E-04                 10  1.34
' 100     195.1      2.40645796828482E-04                11  2.74
' 105     208.8      1.03251578598705E-04                12  6.24
' 110     240.0      4.81695302325482E-05                13 12.0
' 120     360.0                                          14
Select Case P
Case Is > 226.3206     '<11 km
' z = Zs(1) + (Ts(1) / LRs(1)) * (1.-((P / Ps(1)) ^ (-(cRs * LRs(1))/(1000.*cg))))
  z = (44.3307692307692) * (1# - (P / 1013.25) ^ (0.190263235151657))
Case Is > 54.74888     '<20 km
' z = Zs(2) - (cRs * Ts(2) / (cg * 1000)) * ln(P / Ps(2))
  z = 11# - 6.34161998393947 * Log(P / 226.3206)
Case Is > 8.680185     '<32 km
' z = Zs(3) - (Ts(3) / LRs(3)) * (1.-((P / Ps(3)) ^ (-(cRs * LRs(3))/(1000.*cg))))
  z = 20 - 216.65 * (1# - (P / 54.74888) ^ (-2.92712669464088E-02))
Case Is > 1.109063     '<47 km
' z = Zs(4) - (Ts(4) / LRs(4)) * (1.-((P / Ps(4)) ^ ((cRs * LRs(4))/(1000.*cg))))
  z = 32 - 81.6607142857143 * (1# - (P / 8.680185) ^ (-8.19595474499447E-02))
Case Is > 0.6693885    '<51 km
' z = Zs(5) - (cRs * Ts(5) / (cg * 1000)) * ln(P / Ps(5))
  z = 47# - 7.92226839904554 * Log(P / 1.109063)
Case Is > 0.03956419   '<71 km
' z = Zs(6) - (Ts(6) / LRs(6)) * (1.-((P / Ps(6)) ^ ((cRs * LRs(6))/(1000.*cg))))
  z = 51 + 96.6607142857143 * (1# - (P / 0.6693885) ^ (8.19595474499447E-02))
Case Is > 0.003733834  '<84.852 km
' z = Zs(7) - (Ts(7) / LRs(7)) * (1.-((P / Ps(7)) ^ (-(cRs * LRs(7))/(1000.*cg))))
  z = 71 + 107.325 * (1# - (P / 0.03956419) ^ (5.85425338928176E-02))
Case Is > 1.45742511874549E-03   '<90 km
' z = Zs(8) - (cRs * Ts(8) / (cg * 1000)) * ln(P / Ps(8))
  z = 84.852 - 5.47214624555127 * Log(P / 0.003733834)
Case Is > 5.8654139565495E-04       '<95 km
' z = Zs(9) - (Ts(9) / LRs(9)) * (1.-((P / Ps(9)) ^ ((cRs * LRs(9))/(1000.*cg))))
  z = 90# - 642.8679 * (1# - (P / 1.45742511874549E-03) ^ (-8.51208458015383E-03))
Case Is > 2.40645796828482E-04  '<100 km
' z = Zs(4) - (Ts(10) / LRs(10)) * (1.-((P / Ps(10)) ^ ((cRs * LRs(10))/(1000.*cg))))
  z = 95# - 140.597 * (1# - (P / 5.8654139565495E-04) ^ (-3.92234986852218E-02))
Case Is > 1.03251578598705E-04 '<105 km
' z = Zs(4) - (Ts(4) / LRs(4)) * (1.-((P / Ps(4)) ^ ((cRs * LRs(4))/(1000.*cg))))
  z = 100# - 71.20438 * (1# - (P / 2.40645796828482E-04) ^ (-8.02032717123127E-02))
Case Is > 4.81695302325482E-05 '<110 km
' z = Zs(12) - (Ts(12) / LRs(12)) * (1.-((P / Ps(12)) ^ ((cRs * LRs(12))/(1000.*cg))))
  z = 105# - 33.46154 * (1# - (P / 1.03251578598705E-04) ^ (-0.18265269904593))
Case Else
' z = Zs(4) - (Ts(4) / LRs(4)) * (1.-((P / Ps(4)) ^ ((cRs * LRs(4))/(1000.*cg))))
  If P <= 0 Then P = 0.000001
  z = 110 - 20# * (1# - (P / 4.81695302325482E-05) ^ (-0.351255203356906))
End Select
fPtoZ = z
End Function

US Standard Atmosphere Temperature Profile

Function fTstd!(z!)
' Temperature structure of 1976 US Standard Atmosphere
' z in km, fTstd in K
' MJ Mahoney JPL 19980510
Select Case z
Case Is <= 11: fTstd = 288.15 - 6.5 * z
Case Is <= 20: fTstd = 216.65
Case Is <= 32: fTstd = 216.65 + (z - 20)
Case Is <= 47: fTstd = 228.65 + 2.8 * (z - 32)
Case Is <= 51: fTstd = 270.65
Case Is <= 71: fTstd = 270.65 - 2.8 * (z - 51)
Case Is <= 84.852: fTstd = 214.65 - 2# * (z - 71)
Case Is <= 90: fTstd = 186.946
Case Is <= 95: fTstd = 186.946 + 0.2908 * (z - 90)
Case Is <= 100: fTstd = 188.4 + 1.34 * (z - 95)
Case Is <= 105: fTstd = 195.1 + 2.74 * (z - 100)
Case Is <= 110: fTstd = 208.8 + 6.24 * (z - 105)
Case Is <= 120: fTstd = 240 + 12 * (z - 110)
Case Else: fTstd = 360 + 12 * (z - 120)
End Select
End Function