Atmospheric Model
I have to start somewhere, so I may as well start here.
Aerodynamic drag plays a role in determining a rocket’s flight trajectory. In order to calculate drag, I must first calculate the atmospheric density. Atmospheric sciences are a fascinating yet complex domain. Mercifully, for amateur rocketry I do not need a complex atmospheric model to get accurate results. In relative terms, aerodynamic forces involved in rocketry are low compared to forces generated from thrust. As the rocket ascends and atmospheric density exponentially decreases, drag becomes yet less important still.
I have implemented an atmospheric model into traJulia that is based on the US Standard Atmosphere. This is computationally efficient with sufficient accuracy.
Regan’s Hybrid US76–US62 Atmosphere
The 1976 US Standard Atmosphere (US76) is commonly used by aerospace engineers to predict atmospheric conditions at various altitudes. This model is important when considering aerodynamics and flight dynamics. Despite the availability of newer and higher-fidelity methods, US76 remains widely used due to its simplicity, and in many cases it is more than accurate enough. I like it, and have used it extensively over the years.
Unlike the International Standard Atmosphere (ISA), which is controlled by ISO, the reports describing US Standard Atmospheres are public-domain and freely available.
In the book Dynamics of Atmospheric Re-Entry by Regan and Anandakrishnan, an interesting hybrid is described. It uses US76 up to an altitude of 86 km, where the model ends, and extends it using the 1962 US Standard Atmosphere (US62) for altitudes up to 700 km.
For traJulia’s atmospheric model, I have taken inspiration from Regan and implemented something with the same idea, but structured differently to Regan’s original True BASIC code from the 1980s.
Overview of the Code
The atmospheric model function is written in Julia and is structured to calculate five key atmospheric properties at a specified altitude:
- Pressure (Pa)
- Density (kg/m³)
- Temperature (K)
- Speed of sound (m/s)
- Mean free path (m)
In practice, I will likely only ever use density and speed of sound for drag calculations. However, the code also returns pressure, temperature, and mean free path. Regan’s model does this, and these additional terms are computationally cheap in Julia, so it makes sense to include them for potential future use.
The model handles altitudes between 0 and 700 km, which covers the entire mesosphere, thermosphere, and a reasonable portion of the exosphere. Whilst sounding rockets might reach greater altitudes, I do not know of any amateur rockets that approach this. By that point, we are well into the vacuum of space. Drag at these altitudes only matters to spacecraft in orbit for years, not rockets passing through in minutes.
Please check GitHub for the latest version of the script, as it appeared at the time of writing this post.
module AtmosphereReganStandard
"""
This module provides atmospheric properties.
The 1976 US Standard Atmosphere (US76) is used up to 86 km; for altitudes above 86 km
the model is based upon the 1962 US Standard Atmosphere (US62).
This approach is described in:
Regan, F. J., & Anandakrishnan, S. M. (1993). "Dynamics of Atmospheric Re-Entry".
Regan notes:
"The model may easily be changed to meet another standard or to accommodate different
views on the value that should be set for the exospheric temperature. The model
specifies that thermal layers be identified and that within each layer the temperature
varies at most linearly with altitude."
### Arguments:
- `altitude`: Altitude in metres at which to calculate atmospheric properties.
- `T₀`: Ground level temperature in Kelvin. Optional.
- `P₀` Ground level pressure in Pa. Optional.
### Returns:
A named tuple containing:
- `pressure`:Atmospheric pressure at the specified altitude (Pa).
- `density`:Atmospheric density at the specified altitude (kg/m³).
- `temperature`:Kinetic temperature at the specified altitude (K).
- `speed_of_sound`: Speed of sound at the specified altitude (m/s).
- `mean_free_path`: Mean free path length at the specified altitude (m).
### Notes:
- If the specified altitude is outside the range of the model (up to 700 km), the function
returns `missing` for all properties.
- The model assumes that within each atmospheric layer, the temperature varies linearly
with altitude or remains constant (isothermal layer).
"""
using Unitful
using Unitful.DefaultSymbols # For unit symbols like m, K, kg, etc.
# Constants
const R_UNIVERSAL = 8_314.4621u"J/(kmol*K)" # Universal gas constant
const G0 = 9.80665u"m/s^2" # Acceleration due to gravity
const NA = 6.0221e26u"1/kmol" # Avogadro's number
const SIGMA = 3.65e-10u"m" # Effective diameter of atmospheric gas molecule
const M0 = 28.9644u"kg/kmol" # Molecular weight of air at sea level
const RP = 6.3781e6u"m" # Planetary radius
const R_SPECIFIC = R_UNIVERSAL / M0 # Specific gas constant for air
# Atmospheric layer data: (Altitude, Temperature, Molecular Weight)
const DATA = [
(11_000u"m", 216.65u"K", 28.9644u"kg/kmol"),
(20_000u"m", 216.65u"K", 28.9644u"kg/kmol"),
(32_000u"m", 228.65u"K", 28.9644u"kg/kmol"),
(47_000u"m", 270.65u"K", 28.9644u"kg/kmol"),
(51_000u"m", 270.65u"K", 28.9644u"kg/kmol"),
(71_000u"m", 214.65u"K", 28.9644u"kg/kmol"),
(86_000u"m", 186.946u"K", 28.952u"kg/kmol"),
(100_000u"m", 210.65u"K", 28.88u"kg/kmol"),
(110_000u"m", 260.65u"K", 28.56u"kg/kmol"),
(120_000u"m", 360.65u"K", 28.07u"kg/kmol"),
(150_000u"m", 960.65u"K", 26.92u"kg/kmol"),
(160_000u"m", 1110.65u"K", 26.66u"kg/kmol"),
(170_000u"m", 1210.65u"K", 26.5u"kg/kmol"),
(190_000u"m", 1350.65u"K", 25.85u"kg/kmol"),
(230_000u"m", 1550.65u"K", 24.69u"kg/kmol"),
(300_000u"m", 1830.65u"K", 22.66u"kg/kmol"),
(400_000u"m", 2160.65u"K", 19.94u"kg/kmol"),
(500_000u"m", 2420.65u"K", 17.94u"kg/kmol"),
(600_000u"m", 2590.65u"K", 16.84u"kg/kmol"),
(700_000u"m", 2700.00u"K", 16.17u"kg/kmol")
]
# Function to calculate atmospheric properties
function atmosphere(altitude; T₀=288.15, P₀=101325.0)
altitude = Float64(altitude) * u"m" # Convert to Float64 and attach units
T₀ = Float64(T₀) * u"K"
P₀ = Float64(P₀) * u"Pa"
# Check if altitude is within model range
if altitude < 0u"m" || altitude > 700_000u"m"
return (
pressure = missing,
density = missing,
temperature = missing,
speed_of_sound = missing,
mean_free_path = missing
)
end
# Initialise arrays with units
num_layers = length(DATA)
z_breakpoints = [0.0u"m"; [d[1] for d in DATA]]
T_breakpoints = [T₀; [d[2] for d in DATA]]
M_breakpoints = [M0; [d[3] for d in DATA]]
# Compute lapse rates
lapse_rates = diff(T_breakpoints) ./ diff(z_breakpoints)
# Initialise pressure and density arrays
P_breakpoints = zeros(num_layers + 1) .* u"Pa"
rho_breakpoints = zeros(num_layers + 1) .* u"kg/m^3"
P_breakpoints[1] = P₀
rho_breakpoints[1] = P₀ / (R_SPECIFIC * T₀)
# Calculate pressure and density at breakpoints
for i in 1:num_layers
if lapse_rates[i] != 0u"K/m"
# Non-isothermal layer
exponent = -G0 / (R_SPECIFIC * lapse_rates[i])
temperature_ratio = T_breakpoints[i + 1] / T_breakpoints[i]
P_breakpoints[i + 1] = P_breakpoints[i] * temperature_ratio^exponent
rho_breakpoints[i + 1] = rho_breakpoints[i] * temperature_ratio^(exponent - 1)
else
# Isothermal layer
exponent = -G0 * (z_breakpoints[i + 1] - z_breakpoints[i]) / (R_SPECIFIC * T_breakpoints[i])
P_breakpoints[i + 1] = P_breakpoints[i] * exp(exponent)
rho_breakpoints[i + 1] = rho_breakpoints[i] * exp(exponent)
end
end
# Find the atmospheric layer for the given altitude
i = findfirst(x -> altitude < x, z_breakpoints[2:end])
if i === nothing
return (
pressure = missing,
density = missing,
temperature = missing,
speed_of_sound = missing,
mean_free_path = missing
)
end
# Interpolate properties at the desired altitude
if lapse_rates[i] != 0u"K/m"
# Non-isothermal layer
temperature = T_breakpoints[i] + lapse_rates[i] * (altitude - z_breakpoints[i])
exponent = -G0 / (R_SPECIFIC * lapse_rates[i])
temperature_ratio = temperature / T_breakpoints[i]
pressure = P_breakpoints[i] * temperature_ratio^exponent
density = rho_breakpoints[i] * temperature_ratio^(exponent - 1)
else
# Isothermal layer
temperature = T_breakpoints[i]
exponent = -G0 * (altitude - z_breakpoints[i]) / (R_SPECIFIC * T_breakpoints[i])
pressure = P_breakpoints[i] * exp(exponent)
density = rho_breakpoints[i] * exp(exponent)
end
# Interpolate molecular weight
M_start = M_breakpoints[i]
M_end = M_breakpoints[i + 1]
z_start = z_breakpoints[i]
z_end = z_breakpoints[i + 1]
molecular_weight = M_start + (M_end - M_start) * (altitude - z_start) / (z_end - z_start)
# Adjusted specific gas constant
R_specific_adjusted = R_UNIVERSAL / molecular_weight
# Speed of sound
gamma = 1.4 # Ratio of specific heats for air
speed_of_sound = sqrt(gamma * R_specific_adjusted * temperature)
# Mean free path
average_molecular_speed = sqrt(8 * R_SPECIFIC * temperature / π)
collision_frequency = sqrt(2) * π * NA * SIGMA^2 * average_molecular_speed
mean_free_path = average_molecular_speed / collision_frequency
# Return results without units for compatibility
return (
pressure = ustrip(pressure),
density = ustrip(density),
temperature = ustrip(temperature),
speed_of_sound = ustrip(speed_of_sound),
mean_free_path = ustrip(mean_free_path)
)
end
export atmosphere
end
Atmospheric Layers
The model consists of multiple atmospheric layers, each defined by altitude breakpoints, temperature, and molecular weight.
For each layer, a lapse rate (the rate of temperature change with altitude) is calculated. The model then computes pressure and density at the boundaries between layers, determining whether to use an isothermal or non-isothermal formulation. If the lapse rate is non-zero, the temperature varies linearly; otherwise, it remains constant.
The data for these layers, from the troposphere up to the exosphere, are hard-coded into the function, consistent with Regan’s example.
Extensibility
One of the advantages of this approach, as noted by Regan, is that it can easily be adapted to different atmospheric standards or updated to reflect changes in exospheric temperature. Whilst I do not anticipate needing anything “better”, I like the modularity this affords. Even an idiot like me can adjust the parameters with relative ease.
Verification and Validation
The important bit.
Validation is largely unnecessary here. We have decades of confidence in these models, and their validity and limitations are well understood and well outside the scope of this post. If that needs arguing, you are the wrong audience.
Verification, however, is important, particularly as I have relied heavily on LLMs to help shape the code. I have found them useful for building simple code and learning a new language. Julia is new to me. It is a strange experience arguing with an AI about why some decisions are idiotic, while also being impressed by code structures I would not have discovered otherwise. This makes me especially cautious of the output.
For reference data, I have used Digital Dutch’s online 1976 Standard Atmosphere calculator, which I have trusted for many years:
https://www.digitaldutch.com/atmoscalc/
By treating this as an independent and accurate US76 implementation, and using tabulated output, I can compare results.
Verification – Density
As Digital Dutch’s implementation is limited to 86 km, comparison is restricted to this range. I cannot spot any obvious issues with density. Pressure and density would be better plotted on a logarithmic scale, particularly when extending to 700 km, but for the US76 region this looks reasonable.
Verification – Pressure
Pressure shows the same behaviour.
Verification – Speed of Sound and Temperature
Speed of sound agrees well, which, as expected, mirrors the temperature behaviour. There is some divergence at the top of the model that I intend to explore further. The linear variation between breakpoints is consistent with how the code is structured.
I believe Regan’s hybrid model was originally written in imperial units and later converted to metric for publication. This results in some constants, such as the universal gas constant and sea-level gravity, being close to but not exactly equal to modern accepted values. I thought I had corrected all of these, but this discrepancy may indicate something I have missed.
Further Checks
I still need to verify the extended atmosphere above 86 km, as well as the mean free path calculation. The discrepancy in the highest layer needs investigation, and I would like to add some unit testing.
For now, I am tired, and going to bed.