A new Relative-Airmass formula

June 2026

Daniel Pereira Valadés

daniel.pereira.valades@gmail.com

Abstract

This article presents a physically grounded method to compute relative airmass as a function of the apparent solar zenith angle. Unlike traditional phenomenological formulations, the proposed approach models the atmosphere as a finite layer of 100 km thickness, combining a realistic vertical density profile derived from the ideal gas law, hydrostatic equilibrium, and standard atmospheric temperature distributions. The airmass is obtained through numerical integration along the optical path of solar radiation, explicitly accounting for Earth curvature and path geometry. Based on the numerical results, a polynomial approximation is developed to provide a practical and computationally efficient expression. The proposed model shows consistent behavior across a wide range of zenith angles and improves the physical interpretability of airmass calculations, particularly near the horizon. The artcile concludes with two simple equations that can be easily implemented into solar resource tools or irradiance models.

Introduction

The usual formulas for the relative airmass are purely phenomenological. In this article, a novel approach is presented, based on a 100km-thickness atmosphere, to calculate the realtive airmass as a function of the apparent solar zenith angle.

In first place, a model of the atmospheric density is required to determine the relative airmass seen by an observer standing on the surface of the arth.

Atmospheric modelling

An expression for the atmospheric density can be approached by the ideal gas law:

\frac{P}{ρ} = R_{g} \cdot T \Rightarrow ρ = \frac{P}{R_{g} \cdot T}

We will also assume that the temperature is exclusively dependent on the height, or the distance to the earth's surface:

ρ (P, h) = \frac{P}{R_{g} \cdot T (h)}

The temperature can be modeled using the ISO Standard Atmosphere (Table 4 on ISO 2533-1975), where the values of temperature are related to the Geopotential altitude:

Geopotential altitude h (km)	Temperature T (K)
-2	301.15
0	288.15
11	216.65
20	216.65
32	228.65
47	270.65
51	270.65
71	214.65
80	196.65

The pressure can be derived with the hydrostatic equation:

First, calculate the mass contained in the blue cylinder:

d m = ρ \cdot A \cdot d h

The force downwards exerted by that differential of mass is:

d \vec{F} = d m \cdot \vec{g} = ρ \cdot A \cdot d h \cdot \vec{g}

Therefore, the expression for calculating the pressure difference:

P + d P = \frac{d \vec{F} \cdot \vec{n}}{A} = P + \frac{ρ \cdot A \cdot d h \cdot \vec{g} \cdot \vec{n}}{A} = P + ρ \cdot d h \cdot \vec{g} \cdot \vec{n}

Knowing that our criteria is that h grows upwards and P grows downwards, our vectors $\vec{g}$ and $\vec{n}$ are oriented in opposite directions. Therefore:

\vec{g} \cdot \vec{n} = - g

Eq 5 and 6 simplify to:

d P = - ρ \cdot g \cdot d h

Substituting Eq. 2 into Eq. 7, and knowing that the acceleration of gravity depends directly on the height:

d P = - \frac{P}{R_{g} \cdot T (h)} \cdot g (h) \cdot d h

Rearranging terms:

\frac{d P}{P} = - \frac{g (h)}{R_{g} \cdot T (h)} \cdot d h

The expression for the acceleration of gravity can be easily calculated applying Newton's gravity equation:

g = \frac{G \cdot M_{Earth}}{{(R_{Earth} + h)}^{2}}

We can integrate both sides of the equality, to get the expression for calculating the pressure at any given altitude, $l$ :

\int_{P (h = 0)}^{P (h = l)} \frac{d p}{p} = \log (P (l)) - \log (P (0)) = \log (\frac{P (l)}{P_{0}}) = - \frac{G \cdot M_{Earth}}{R_{g}} \cdot \int_{0}^{l} \frac{d h}{T (h) \cdot {(R_{Earth} + h)}^{2}}

The right side of the equality will be solved numerically, since there is no expression for the Temperature as a function of height. The table shown above, plus an extra point at 100km with a temperature of 195K [Source: WolframAlpha, 25-Mar-2026], will be converted into a continuous curve with a modified-Akima cubic interpolation:

We get a numerical expression for the pressure at any given altitude $l$ :

P (l) = P_{0} \cdot \exp (- \frac{G \cdot M_{Earth}}{R_{g}} \cdot \int_{0}^{l} \frac{d h}{T (h) \cdot {(R_{Earth} + h)}^{2}})

And using Eq. 2 into Eq. 12 we get the expression for the density at any given altitude $l$ :

ρ (l) = \frac{P_{0}}{R_{g} \cdot T (l)} \cdot \exp (- \frac{G \cdot M_{Earth}}{R_{g}} \cdot \int_{0}^{l} \frac{d h}{T (h) \cdot {(R_{Earth} + h)}^{2}})

This expression on Eq. 13 can be numerically calculated for any given altitude, using the following values for the constants:

$G = 6.674 \cdot 10^{−11} m^{3} \cdot s^{−2} \cdot {kg}^{−1}$
$M_{Earth} = 5.9722 \cdot 10^{24} kg$
$R_{g} = 287.052874 J \cdot {kg}^{−1} \cdot K^{−1}$
$R_{Earth} = 6371 \cdot 10^{3} m$
$P_{0} = 101325 Pa$

Calculating the length of the column of atmosphere

Now that we can calculate the density of air at any given altitude, we must calculate the length of the column of atmosphere that a solar ray must traverse in order to reach an observer standing on the surface of the Earth.

Let us start the conception of the problem with an observer standing on the surface of the Earth, looking directly at the sun (don't ever do this without the adeqcuate protection: the human eye is extremely sensitive and exposure to direct sunlight can lead to ocular injuiries including blindness).

Assume that the Earth is a perfect sphere and that the sun is very low over the horizon, presenting thus a high zenithal angle θ (notice that in the figure below the Earth and atmosphere sizes are not to scale):

Let us then add some measurements, assuming a constant total atmospheric height H (which we will assume to be 100km) and R for the radius of Earth (6371km as stated above). Also, the vector pointing the sun (in parametric form) and the equation for the top of the atmosphere is presented in the following figure. Note that we assume the origin for the coordinates at the Earth's center:

We can calculate the value of t when the vector intercepts the limit of the atmosphere:

\vec{v} = (x, y) = (t \cdot \sin θ, R + t \cdot \cos θ)

Using the equation of the top of the atmosphere and substituting Eq. 14 into it:

x^{2} + y^{2} = {(R + H)}^{2} \Rightarrow {(t_{\max} \cdot \sin θ)}^{2} + {(R + t_{\max} \cdot \cos θ)}^{2} = {(R + H)}^{2}

This gives the value of $t_{\max}$ , which is equal to the total length of teh column of atmosphere traversed by a sunray reaching the surface of the earth:

t_{\max} (θ) = R \cdot (- \cos θ + \sqrt{{(1 + \frac{H}{R})}^{2} - \sin^{2} θ})

Note that Eq. 16 does not result into infinty when plugging 90° for the zenithal angle, compared to the simple secant formula (that assumes earth is flat). Using R=6371 km, H=100 km and θ=90°:

t_{\max} (90°) = 1133.2255 km

Now that we know the maximum extent of the vector, for simplicity, we will assume that the parametric expression of the vector pointing the sun is:

\vec{v} = (x, y) = (t_{\max} \cdot ξ \cdot \sin θ, R + t_{\max} \cdot ξ \cdot \cos θ)

With ξ ranging from 0 (at the Earth's surface) to 1 (at the top of the atmosphere):

The altitude for any point of the vector can be calculated using the cartesian coordinates, taking advantage of the fact that we took our origin at the center of the Earth.

The distance to the center of the Earth can be easily calculated for any point with ξ from 0 to 1:

d (ξ, θ) = \sqrt{x^{2} + y^{2}} = \sqrt{{(t_{max} (θ) \cdot ξ \cdot \sin θ)}^{2} + {(R + t_{max} (θ) \cdot ξ \cdot \cos θ)}^{2}}

So the altitude at which any point stands over the surface of the Earth is:

l (ξ, θ) = d (ξ, θ) - R

Applying Eq. 18 into Eq. 19, we get the expression for the altitude:

l (ξ, θ) = \sqrt{{(t_{max} (θ) \cdot ξ \cdot \sin θ)}^{2} + {(R + t_{max} (θ) \cdot ξ \cdot \cos θ)}^{2}} - R

Analyzing the last equation, it can be seen that negative altitudes are obtained when θ>90°. This is caused by the fact that a zenithal angle greater than 90° means that the sun is below the horizon. Therefore, a sunray reaching us from below the horizon is only possible when the altiude of the surroundings is lower than the observer altitudes.

The formula to compute the density of air does not change if the altitude is negative (at negative altitudes, the atmosphere is thicker and this will result in higher densities). However, the formula for altitude gives values well below 2km at values over 91.44°. Since our density function is based on the temperature table, and this is limited to a minimum altitude of −2km, it is not recommended to go below this minumim altitude. Even more, it would not be recommended to go even below 0, since no sun-rays will go through that direction (but rather be diffracted through low elevations).

Considering this, Eq. 21 should be bounded to the range 0-100km. Nevertheless, to show that both the altitude and the density functions work even with high zenith angles, the following graphs are shown:

Note that the altitude $l$ follows the expected path:

for low zenithal angles it almost grows linearly with the parameter ξ,
for θ>90° it becomes negative until certain point (when it reaches the surface of the sphere),
and every curve starts at 0 (ξ=0 is on the ground) and finishes at 100km (ξ=1 is just where the atmosphere ends, whcih we defined as 100km).

Similarly, the density function follows the expected path, given the altitudes for each path.

The corrected equation, with a lower limit of 0km altitude is, therefore:

l (ξ, θ) = max (0, \sqrt{{(t_{max} (θ) \cdot ξ \cdot \sin θ)}^{2} + {(R + t_{max} (θ) \cdot ξ \cdot \cos θ)}^{2}} - R)

The resulting graphs are:

With altitudes ≥0 and densities limited to ρ₀ (density at 0m altitude).

Combination of atmospheric modelling and length of the column

Now, we have a parameter $ξ$ , which we know must change from 0 to 1. Also, we have the altitude relative to the ground $l (ξ)$ [Eq. 22]. We have the density of the air as a function of the altitude $ρ (l)$ [Eq. 13].

To do the numeric integration, we can imageine that we are a particle with area A, travelling from ξ=0 to ξ=1, calculating at each step the altitude at our current position inside the path, the density of the air at this position and the distance differential (dx) we have traveled from the previous value of ξ. This gives us a volume (our area A times the distance differential traveled) to be multiplied by the density to calculate the mass in that portion. The distance differential is the total distance that we must travel (t_max) multiplied by the differential of the variable that goes from 0 to 1 (dξ). Summing all the portions will give us the total air mass. Calculating everything regardless the value of A, the result will be a mass by unit area (kg/m²):

dm = ρ (l (ξ, θ)) \cdot A \cdot dx = ρ (l (ξ, θ)) \cdot A \cdot t_{max} (θ) \cdot d ξ

Summing all the steps, and using a very small step size (1·10^-5) allows us to obtain the total mass per unit area:

\frac{m (θ)}{A} = \frac{1}{A} \cdot \sum_{ξ = 0}^{1} dm = t_{max} (θ) \cdot \sum_{ξ = 0}^{1} ρ (l (ξ, θ)) \cdot d ξ

The relative airmass will be:

AM (θ) = \frac{m (θ)}{m (0)}

Results

Using Eq. 24 and 25 with the formula for altitude l [Eq. 22], the formula for density ρ [Eq. 13], and the formula for the total column height t_max [Eq. 16] we get the following results for AM(θ), by performing the abovementioned numeric integration with step 10^-5 with zenithal angles θ ranging from 0° to 96°.

Approximation formula

Since the model is complicated to calculate, an approximation formula is going to be derived for the results shown above. The approximation to these results should be such that the formula is easy to implement and that the error between the data obtained and the approximant is neglibile.

An approach for this calculation is checking the most influential parameter on the value of airmass.

Having studied how the parameters affect, it is very obvious that the most relevant is the length of the column of atmosphere, t_max. In fact, having considered a constant atmospheric thickness, creates a direct correlation between the total column height and the path of a sunray reaching the surface of the earth.

The influence of the value of t_max over the relative airmass AM is shown below:

It can be seen that there is a very strong correlation between the two variables. A polynomial approximant can be used to calculate this relationship (note that the graph has logarithmic scale in the Y axis). The coefficients for the polynomial are obtained by least squares method, with a polynomial of degree 7 correlating t_max and log(AM). In order to obtain coefficients in the same order of magnitude, the polynomial will be calculated with $\frac{t_{\max}}{R}$ .

X = [\begin{matrix} 1 & \frac{{t_{\max}}_{1}}{R} & {(\frac{{t_{\max}}_{1}}{R})}^{2} & {(\frac{{t_{\max}}_{1}}{R})}^{3} & {(\frac{{t_{\max}}_{1}}{R})}^{4} & {(\frac{{t_{\max}}_{1}}{R})}^{5} & {(\frac{{t_{\max}}_{1}}{R})}^{6} & {(\frac{{t_{\max}}_{1}}{R})}^{7} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & \frac{{t_{\max}}_{i}}{R} & {(\frac{{t_{\max}}_{i}}{R})}^{2} & {(\frac{{t_{\max}}_{i}}{R})}^{3} & {(\frac{{t_{\max}}_{i}}{R})}^{4} & {(\frac{{t_{\max}}_{i}}{R})}^{5} & {(\frac{{t_{\max}}_{i}}{R})}^{6} & {(\frac{{t_{\max}}_{i}}{R})}^{7} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & \frac{{t_{\max}}_{n}}{R} & {(\frac{{t_{\max}}_{n}}{R})}^{2} & {(\frac{{t_{\max}}_{n}}{R})}^{3} & {(\frac{{t_{\max}}_{n}}{R})}^{4} & {(\frac{{t_{\max}}_{n}}{R})}^{5} & {(\frac{{t_{\max}}_{n}}{R})}^{6} & {(\frac{{t_{\max}}_{n}}{R})}^{7} \end{matrix}]

Y = [\begin{matrix} log ({AM}_{1}) \\ ⋮ \\ log ({AM}_{i}) \\ ⋮ \\ log ({AM}_{n}) \end{matrix}]

\hat{β} = {(X' \cdot X)}^{-1} \cdot {(X' \cdot Y)}^{}

The value obtained for $\hat{β}$ for the data shown above is:

\hat{β} = [\begin{matrix} -1.165739 \\ 96.068758 \\ -1598.662757 \\ 17276.030762 \\ -109590.458437 \\ 414954.803058 \\ -867699.500142 \\ 762235.787430 \end{matrix}]

Resulting into the expression:

\log (AM) = −1.165739 +96.068758 \cdot \frac{t_{\max}}{R} −1598.662757 \cdot {(\frac{t_{\max}}{R})}^{2} +17276.030762 \cdot {(\frac{t_{\max}}{R})}^{3} −109590.458437 \cdot {(\frac{t_{\max}}{R})}^{4} +414954.803058 \cdot {(\frac{t_{\max}}{R})}^{5} −867699.500142 \cdot {(\frac{t_{\max}}{R})}^{6} +762235.787430 \cdot {(\frac{t_{\max}}{R})}^{7}

The value of $\frac{t \max}{R}$ can be simplified, using Eq.16:

\frac{t_{\max} (θ)}{R} = - \cos θ + \sqrt{{(1 + \frac{H}{R})}^{2} - \sin^{2} θ}

This can be simplified to:

\frac{t_{\max} (θ)}{R} = - \cos θ + \sqrt{{(1 + \frac{H}{R})}^{2} - 1 + \cos^{2} θ}

Using R = 6371 km and H = 100km:

\frac{t_{\max} (θ)}{R} = - \cos θ + \sqrt{0.0316386144 + \cos^{2} θ}

Then using Eq. 30 with this value of $\frac{t \max}{R}$ , we can compare to the curve for AM obtained by numeric integration:

Final expression and Comparison to other models

The model proposed in this article:

$AM = exp (−1.165739 +96.068758 \cdot τ −1598.662757 \cdot τ^{2} +17276.030762 \cdot τ^{3} −109590.458437 \cdot τ^{4} +414954.803058 \cdot τ^{5} −867699.500142 \cdot τ^{6} +762235.787430 \cdot τ^{7})$
Where $τ = \frac{t_{\max}}{R} = - \cos θ + \sqrt{0.0316386144 + \cos^{2} θ}$

will be compared against other known airmass models:

sec(θ):

$AM = \frac{1}{\cos θ}$

Kasten, 1966:

$AM = \frac{1}{cos θ +0.15⋅ {(93.885− θ)}^{−1.253}}$

Young-Irvine, 1967:

$AM = \frac{1}{cos θ ⋅ (1 − 0.0012 ⋅ (\frac{1}{\cos^{2} θ} − 1))}$

Kasten-Young, 1989:

$AM = \frac{1}{cos θ +0.50572⋅ {(96.07995− θ)}^{−1.6364}}$

Gueymard, 1993:

$AM = \frac{1}{cos θ + 0.00176759 ⋅ θ ⋅ {(94.37515− θ)}^{−1.21563}}$

Young, 1994:

$AM = \frac{1.002432 ⋅ {cos}^{2} θ + 0.148386 ⋅ cos θ + 0.0096467}{{cos}^{3} θ + 0.149864 ⋅ {cos}^{2} θ + 0.0102963 ⋅ cos θ + 0.000303978}$

Pickering, 2002:

$AM = \frac{1}{sin (90 − θ + \frac{244}{165+47⋅ {(90- θ)}^{1.1}})}$

Results for θ ∈ [0,96°]:

Zooming near θ = 90°:

Conclusion

The model developed in this work provides a physically consistent framework for the calculation of relative airmass by explicitly coupling atmospheric density profiles with the geometrical path of solar radiation through a finite atmosphere. By combining a realistic atmosphere description based on thermodynamic principles with a geometric formulation that accounts for Earth curvature, the resulting airmass values remain well-defined even at high zenith angles, avoiding the divergences associated with simpler formulations such as the secant approximation.

The numerical integration approach allows the airmass to be interpreted as the total mass of air traversed per unit area, offering a clear physical meaning that is often lacking in empirical models. Furthermore, the strong correlation identified between the airmass and the normalized atmospheric path length enables the derivation of an accurate polynomial approximation, significantly simplifying practical implementation while preserving high fidelity with respect to the full model.

Comparisons with established models show that the proposed formulation captures the expected behavior across the full angular range, with particularly improved performance near the horizon, where many classical approximations lose validity. This makes the model especially relevant for applications where low solar elevations are critical, such as solar energy resource assessment and atmospheric optics.

Future work may focus on extending the model to account for observer altitude, allowing for more accurate representation of elevated sites and enabling additional applications such as horizon distance estimation. Incorporating more detailed atmospheric profiles or temporal variability could further enhance the model’s applicability in high-precision studies.

Further developments

A potential improvement regarding the observer actual altitude can be made, by including the observer altitude on Eq. 14. The vector in parametric form could have been (t⋅sinθ,R+Alt+t⋅cosθ). This also affects equations 15 to 19. Inclduing this Alt value will also allow to calculate the distance to the horizon at a given altitude, by equating l(ξ)=0 and finding the value of ξ⋅t_max. However, and despite improving the model, this would mean that the approximated expression for the relative airmass will be different for every altitude (or would need extra adjustment to account for the altitude).

The actual shape of the Earth using the reference geoid can also be considered in further works, instead of assuming a spherical Earth, therefore including the actual location in the calculation of the relative airmass.

Content