Repository: cta-observatory/ctapipe
Pull request #2000 (2022-07-08 15:44:57)

Atmosphere models for h_max to X_max conversion
### Features

Adds an `AtmosphereDensityProfile` class hierarchy with three implementations:
- `ExponentialAtmosphereDensityProfile`  (purely mathematical)
- `TableAtmsophereDensityProfile` (cubic interpolation of a tabulated model)
- `FiveLayerAtmosphereDensityProfile`  (implements the Corsika 5-layer representation where each layer is either an exponential or linear model)

This adds a new `atmosphere_density_profile`  attribute to EventSource, with an implementation in `SimTelEventSource` (reading them from simtel files) and `HDF5EventSource`, as well as in `DataWriter` (propagating them to output files).

By default SimTelEventSource will return a TableAtmosphereDensityProfile if it is found in the data, otherwise it falls back to a FiveLayerAtmosphereDensityProfile, or None if no profile exists.  This is user-selectable with the `SimTelEventSource.atmosphere_model_choice` option (AUTO, TABLE, FIVELAYER)

### Example:

```py3
from ctapipe.atmosphere import ExponentialAtmosphereDensityProfile, TableAtmosphereDensityProfile

# basic profile, using default parameters
print(density)
print("Density:       ", density(10*u.km))
print("Column Density (0d):", density.line_of_sight_integral(distance=10*u.km, zenith_angle=0*u.deg))
print("Column Density (45d):", density.line_of_sight_integral(distance=10*u.km, zenith_angle=45*u.deg))
```
```
ExponentialAtmosphereDensityProfile(h0=<Quantity 8. km>, rho0=<Quantity 0.00125 g / cm3>)
Density:        0.0003581309960752376 g / cm3
Column Density (0d): 286.50479686019 g / cm2
Column Density (45d): 584.3180219705157 g / cm2
```

Loading one from a simtel or ctapipe file can be done with an EventSource:
```py3
# Table profile loaded from a SimTel file
with EventSource("gammas.simtel.zst", atmosphere_profile_choice="TABLE") as source:
     density_tab = source.atmosphere_density_profile
print(density_tab)
```

```
TableAtmosphereDensityProfile(meta={'obs_level': <Quantity 215800. cm>, 'atmo_id': 99, 
    'atmo_name': b'atmprof_ecmwf_north_winter_fixed.dat', 'htoa': 12000000.0}, rows=55)
```

All profiles can also be plotted:
```py3
density.peek()
density_tab.peek()
```
<img width="1036" alt="image" src="https://user-images.githubusercontent.com/11677812/178535922-0590a560-bc18-40de-84c7-cd11c38e4c6a.png">





review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/cta-observatory/ctapipe/pull/2000"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>
codecov[bot]: # [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
> Merging [#2000](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) (329a1e4) into [master](https://codecov.io/gh/cta-observatory/ctapipe/commit/2b08140df8ac623c6bde09402dd918a7999db830?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) (2b08140) will **decrease** coverage by `0.15%`.
> The diff coverage is `68.26%`.

> :exclamation: Current head 329a1e4 differs from pull request most recent head ddaace4. Consider uploading reports for the commit ddaace4 to get more accurate results

```diff
@@            Coverage Diff             @@
##           master    #2000      +/-   ##
==========================================
- Coverage   92.35%   92.19%   -0.16%     
==========================================
  Files         193      196       +3     
  Lines       15585    15689     +104     
==========================================
+ Hits        14393    14464      +71     
- Misses       1192     1225      +33     
```


| [Impacted Files](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | Coverage Δ | |
|---|---|---|
| [ctapipe/atmosphere/model.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9hdG1vc3BoZXJlL21vZGVsLnB5) | `61.17% <61.17%> (ø)` | |
| [ctapipe/atmosphere/\_\_init\_\_.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9hdG1vc3BoZXJlL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [ctapipe/atmosphere/tests/test\_models.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9hdG1vc3BoZXJlL3Rlc3RzL3Rlc3RfbW9kZWxzLnB5) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Last update [2b08140...ddaace4](https://codecov.io/gh/cta-observatory/ctapipe/pull/2000?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).

kbernloehr: Hi Max,

I am not aware of simtel files from a single run with more than
an I/O block of type 1216 - if you found such a case, let me
know how it was produced. If they are identical, using the
first one is OK but there are situations possible where a
first one might just contain part of the information and
might later get complemented by a second such block
with the other part. Typical CORSIKA simulation should
fill both parts right away though.

So, well, there are two parts of content in this data block.
One part, which you probably want to convert, is the copy
of the atmprof*.dat file read from CORSIKA. The other part
is the result of the 5-layer fit used by the CORSIKA EGS part.
The code running with CORSIKA should report both parts,
if available, but you might end up with just one part.
The table part would be missing if you run corsika without
an atmospheric profile table. The other ways to configure
the atmospheric profiles in CORSIKA only set the 5-layer
parameters. In that case there is no table to report.
The thickc() and heighc() functions in mc_atmprof.c are
C implementations of what CORSIKA is doing internally
with the 5-layer atmosphere models.

Too bad that EGS part is hard-coded to use the 5-layer
function but it cannot accurately describe the actual
atmosphere. If you want to lookup some height from an
atmospheric depth in the longitudinal profiles, you first
need to make up your mind which description is the more
appropriate one for the task. Since most Cherenkov light
gets emitted from e+-, the 5-layer description might actually
be closer to what CORSIKA originally generated. Xmax is a
mixed beast, at least for hadron showers, so there would be
no really correct way. Hadrons and muons follow directly
the table profile, if any got set up and you don't run a
CURVED option CORSIKA.

Note that interpolation of the tabulated profiles is quite non-trivial.
In particular, you want to have height(thick(h)) == h, at an
acceptable accuracy.
Linear interpolation in a mostly exponential density profile
will obviously give wrong answers, even though that looks
reversible. The atmo.c code goes through some effort to
cspline interpolate in log(thick) versus height - except at the
top where thick=0 is not suitable for the log while csplines
are not very forgiving for any sudden transitions (jumps in
second derivative). The reversion through the heighx_ and
thickx_ functions in atmo.c are good to about 0.2 mm throughout
the whole atmosphere while direct interpolation methods
might easily miss by a 100 meters.

   Best regards,
     Konrad

Am 11.07.22 um 18:08 schrieb Maximilian Nöthe:
> @maxnoe commented on this pull request.
>
>
>
>> +
> +
> +def read_simtel_profile(simtelfile: str) -> Table:
> +    """
> +    Read an atmosphere profile from a SimTelArray file as an astropy Table
> +
> +    Returns
> +    -------
> +    Table:
> +        table with columns `height`, `density`, and `column_density`
> +        along with associated metadata
> +    """
> +    import eventio
> +
> +    tables = []
> +    with eventio.SimTelFile(simtelfile) as simtel:
>
> @kbernloehr Any reason why there could be more than one atmospheric profile in a simtel file?
>


kosack: Thanks.  As it's not intended for simulation, the use case for this package probably doesn't require extreme precision in the atmosphere profile or it's interpolation as other systematics dominate, but in any case I'll switch to a log  interpolation and use splines.  I'll also add a disclaimer in the documentation!  
kbernloehr: Hi Karl,

I agree that for Xmax lookup you can be a bit more
relaxed. You definitely want to leave out the last point
(thick=0) but perhaps even the second to last point
to avoid wiggles in the csplines. The shower maximum
will never be at these altitudes.

And, even though all of our available productions used
tabulated atmospheres, you might want to be prepared
for those where the atmospheric profile is defined by other
means (switch to 5-layer description).

    Best,
     Konrad

Am 12.07.22 um 10:11 schrieb Karl Kosack:
> Thanks.  As it's not intended for simulation, the use case for this package probably doesn't require extreme precision in the atmosphere profile or it's interpolation as other systematics dominate, but in any case I'll switch to a log  interpolation and use splines.  I'll also add a disclaimer in the documentation!
>


maxnoe: > means (switch to 5-layer description).

Do you want to add a class for that description as well, @kosack ? As a fallback if the tabulated version is not available?
kosack: Interpolation in the wiggly areas looks like this now (with cubic splines used):

<img width="618" alt="image" src="https://user-images.githubusercontent.com/11677812/178747376-bf52e4a4-9cfd-4005-9fe8-20f6c04c7483.png">

kosack: >> means (switch to 5-layer description).

> Do you want to add a class for that description as well, @kosack ? As a fallback if the tabulated version is not available?

Is that documented somewhere and do we have a file that has that format atmosphere? If so, I can add it, or  we could always add a subclass in a later version
kosack: I added a `FiveLayer` profile. Not sure I did it correctly, but at least the column density profile looks ok.
kosack: The last item to implement here is the i/o.  Since there could in principle be more than one version of the atmosphere profile in the file (one from the simulation, and one from real measurements for observation data), I think the ones loaded by the SimTelEventSource will go in */simulation/service/atmosphere_profile*, which then allows one to add a */dl2/service/atmosphere_profile* later for the frequently updated/measured one.
kosack: Strange, I get nearly, but not exactly the same result as @kbernloehr  using the data above:

Thickness is identical:
<img width="568" alt="image" src="https://user-images.githubusercontent.com/11677812/180001071-8979679c-f159-4447-997a-12771820bdc3.png">

Rho is nearly identical, except at one transition 
<img width="557" alt="image" src="https://user-images.githubusercontent.com/11677812/180000949-7d21496e-53d1-41fc-b287-5f83f8c03e48.png">


Patch:
diff --git a/ctapipe/atmosphere.py b/ctapipe/atmosphere.py
new file mode 100644
index 00000000..f544c634
--- /dev/null
+++ b/ctapipe/atmosphere.py
@@ -0,0 +1,396 @@
+#!/usr/bin/env python3
+
+"""Atmosphere density models and functions to transform between column density
+(X in grammage units) and height (meters) units.
+
+Zenith angle is taken into account in the line-of-sight integral to compute the
+column density X assuming Earth as a flat plane (the curvature is not taken into
+account)
+
+"""
+
+import abc
+from dataclasses import dataclass
+from functools import partial
+from typing import Dict
+
+import numpy as np
+from astropy import units as u
+from astropy.table import Table
+from scipy.interpolate import interp1d
+
+__all__ = [
+    "AtmosphereDensityProfile",
+    "ExponentialAtmosphereDensityProfile",
+    "TableAtmosphereDensityProfile",
+    "FiveLayerAtmosphereDensityProfile",
+]
+
+SUPPORTED_TABLE_VERSIONS = {
+    1,
+}
+
+
+class AtmosphereDensityProfile(abc.ABC):
+    """
+    Base class for models of atmosphere density.
+    """
+
+    @abc.abstractmethod
+    def __call__(self, height: u.Quantity) -> u.Quantity:
+        """
+        Returns
+        -------
+        u.Quantity["g cm-3"]
+            the density at height h
+        """
+
+    @abc.abstractmethod
+    def integral(self, height: u.Quantity) -> u.Quantity:
+        r"""Integral of the profile along the height axis, i.e. the *atmospheric
+        depth* :math:`X`.
+
+        .. math:: X(h) = \int_{h}^{\infty} \rho(h') dh'
+
+        Returns
+        -------
+        u.Quantity["g/cm2"]:
+            Integral of the density from height h to infinity
+
+        """
+
+    def line_of_sight_integral(
+        self, distance: u.Quantity, zenith_angle=0 * u.deg, output_units=u.g / u.cm**2
+    ):
+        r"""Line-of-sight integral from the shower distance to infinity, along
+        the direction specified by the zenith angle. This is sometimes called
+        the *slant depth*. The atmosphere here is assumed to be Cartesian, the
+        curvature of the Earth is not taken into account.
+
+        .. math:: X(h, \Psi) = \int_{h}^{\infty} \rho(h' \cos{\Psi}) dh'
+
+        Parameters
+        ----------
+        distance: u.Quantity["length"]
+           line-of-site distance from observer to point
+        zenith_angle: u.Quantity["angle"]
+           zenith angle of observation
+        output_units: u.Unit
+           unit to output (must be convertible to g/cm2)
+
+        """
+
+        return (
+            self.integral(distance * np.cos(zenith_angle)) / np.cos(zenith_angle)
+        ).to(output_units)
+
+    def peek(self):
+        """
+        Draw quick plot of profile
+        """
+        # pylint: disable=import-outside-toplevel
+        import matplotlib.pyplot as plt
+
+        fig, axis = plt.subplots(1, 3, constrained_layout=True, figsize=(10, 3))
+
+        fig.suptitle(self.__class__.__name__)
+        height = np.linspace(1, 100, 500) * u.km
+        density = self(height)
+        axis[0].set_xscale("linear")
+        axis[0].set_yscale("log")
+        axis[0].plot(height, density)
+        axis[0].set_xlabel(f"Height / {height.unit.to_string('latex')}")
+        axis[0].set_ylabel(f"Density / {density.unit.to_string('latex')}")
+        axis[0].grid(True)
+
+        distance = np.linspace(1, 100, 500) * u.km
+        for zenith_angle in [0, 40, 50, 70] * u.deg:
+            column_density = self.line_of_sight_integral(distance, zenith_angle)
+            axis[1].plot(distance, column_density, label=f"$\\Psi$={zenith_angle}")
+
+        axis[1].legend(loc="best")
+        axis[1].set_xlabel(f"Distance / {distance.unit.to_string('latex')}")
+        axis[1].set_ylabel(f"Column Density / {column_density.unit.to_string('latex')}")
+        axis[1].set_yscale("log")
+        axis[1].grid(True)
+
+        zenith_angle = np.linspace(0, 80, 20) * u.deg
+        for distance in [0, 5, 10, 20] * u.km:
+            column_density = self.line_of_sight_integral(distance, zenith_angle)
+            axis[2].plot(zenith_angle, column_density, label=f"Height={distance}")
+
+        axis[2].legend(loc="best")
+        axis[2].set_xlabel(
+            f"Zenith Angle $\\Psi$ / {zenith_angle.unit.to_string('latex')}"
+        )
+        axis[2].set_ylabel(f"Column Density / {column_density.unit.to_string('latex')}")
+        axis[2].set_yscale("log")
+        axis[2].grid(True)
+
+        plt.show()
+        return fig, axis
+
+    @classmethod
+    def from_table(cls, table: Table):
+        """return a subclass of AtmosphereDensityProfile from a serialized
+        table"""
+
+        if "TAB_TYPE" not in table.meta:
+            raise ValueError("expected a TAB_TYPE metadata field")
+
+        version = table.meta.get("TAB_VER", "")
+        if version not in SUPPORTED_TABLE_VERSIONS:
+            raise ValueError(f"Table version not supported: '{version}'")
+
+        tabtype = table.meta.get("TAB_TYPE")
+
+        if tabtype == "ctapipe.atmosphere.TableAtmosphereDensityProfile":
+            return TableAtmosphereDensityProfile(table)
+        if tabtype == "ctapipe.atmosphere.FiveLayerAtmosphereDensityProfile":
+            return FiveLayerAtmosphereDensityProfile(table)
+
+        raise TypeError(f"Unknown AtmosphereDensityProfile type: '{tabtype}'")
+
+
+@dataclass
+class ExponentialAtmosphereDensityProfile(AtmosphereDensityProfile):
+    """
+    A simple functional density profile modeled as an exponential.
+
+    The is defined following the form:
+
+    .. math:: \\rho(h) = \\rho_0 e^{-h/h_0}
+
+
+    .. code-block:: python
+
+        from ctapipe.atmosphere import ExponentialAtmosphereDensityProfile
+        density_profile = ExponentialAtmosphereDensityProfile()
+        density_profile.peek()
+
+
+    Attributes
+    ----------
+    scale_height: u.Quantity["m"]
+        scale height (h0)
+    scale_density: u.Quantity["g cm-3"]
+        scale density (rho0)
+    """
+
+    scale_height: u.Quantity = 8 * u.km
+    scale_density: u.Quantity = 0.00125 * u.g / (u.cm**3)
+
+    @u.quantity_input(height=u.m)
+    def __call__(self, height) -> u.Quantity:
+        return self.scale_density * np.exp(-height / self.scale_height)
+
+    @u.quantity_input(height=u.m)
+    def integral(
+        self,
+        height,
+    ) -> u.Quantity:
+        return (
+            self.scale_density * self.scale_height * np.exp(-height / self.scale_height)
+        )
+
+
+class TableAtmosphereDensityProfile(AtmosphereDensityProfile):
+    """Tabular profile from a table that has both the density and it's integral
+    pre-computed.  The table is interpolated to return the density and its integral.
+
+    .. code-block:: python
+
+        from astropy.table import Table
+        from astropy import units as u
+
+        from ctapipe.atmosphere import TableAtmosphereDensityProfile
+
+        table = Table(
+            dict(
+                height=[1,10,20] * u.km,
+                density=[0.00099,0.00042, 0.00009] * u.g / u.cm**3
+                column_density=[1044.0, 284.0, 57.0] * u.g / u.cm**2
+            )
+        )
+
+        profile = TableAtmosphereDensityProfile(table=table)
+        print(profile(10 * u.km))
+
+
+    Attributes
+    ----------
+    table: Table
+        Points that define the model
+
+    See Also
+    --------
+    ctapipe.io.eventsource.EventSource.atmosphere_density_profile:
+        load a TableAtmosphereDensityProfile from a supported EventSource
+    """
+
+    def __init__(self, table: Table):
+        """
+        Parameters
+        ----------
+        table: Table
+            Table with columns `height`, `density`, and `column_density`
+        """
+
+        for col in ["height", "density", "column_density"]:
+            if col not in table.colnames:
+                raise ValueError(f"Missing expected column: {col} in table")
+
+        self.table = table[
+            (table["height"] >= 0)
+            & (table["density"] > 0)
+            & (table["column_density"] > 0)
+        ]
+
+        # interpolation is done in log-y to minimize spline wobble
+
+        self._density_interp = interp1d(
+            self.table["height"].to("km").value,
+            np.log10(self.table["density"].to("g cm-3").value),
+            kind="cubic",
+        )
+        self._col_density_interp = interp1d(
+            self.table["height"].to("km").value,
+            np.log10(self.table["column_density"].to("g cm-2").value),
+            kind="cubic",
+        )
+
+        # ensure it can be read back
+        self.table.meta["TAB_TYPE"] = "ctapipe.atmosphere.TableAtmosphereDensityProfile"
+        self.table.meta["TAB_VER"] = 1
+
+    @u.quantity_input(height=u.m)
+    def __call__(self, height) -> u.Quantity:
+        return u.Quantity(
+            10 ** self._density_interp(height.to_value(u.km)), u.g / u.cm**3
+        )
+
+    @u.quantity_input(height=u.m)
+    def integral(self, height) -> u.Quantity:
+        return u.Quantity(
+            10 ** self._col_density_interp(height.to_value(u.km)), u.g / u.cm**2
+        )
+
+    def __repr__(self):
+        return (
+            f"{self.__class__.__name__}(meta={self.table.meta}, rows={len(self.table)})"
+        )
+
+
+# Here we define some utility functions needed to build the piece-wise 5-layer
+# model.
+
+# pylint: disable=invalid-name,unused-argument
+def _exponential(h, a, b, c):
+    """exponential atmosphere"""
+    return a + b * np.exp(-h / c)
+
+
+def _d_exponential(h, a, b, c):
+    """derivative of exponential atmosphere"""
+    return -b / c * np.exp(-h / c)
+
+
+def _linear(h, a, b, c):
+    """linear atmosphere"""
+    return a - b * h / c
+
+
+def _d_linear(h, a, b, c):
+    """derivative of linear atmosphere"""
+    return -b / c
+
+
+class FiveLayerAtmosphereDensityProfile(AtmosphereDensityProfile):
+    r"""
+    CORSIKA 5-layer atmosphere model
+
+    Layers 1-4  are modeled with:
+
+    .. math:: T(h) = a_i + b_i \exp{-h/c_i}
+
+    Layer 5 is modeled with:
+
+    ..math:: T(h) = a_5 - b_5 \frac{h}{c_5}
+
+    References
+    ----------
+    [corsika-user] D. Heck and T. Pierog, "Extensive Air Shower Simulation with CORSIKA:
+        A User’s Guide", 2021, Appendix F
+    """
+
+    def __init__(self, table: Table):
+        self.table = table
+
+        param_a = self.table["a"].to("g/cm2")
+        param_b = self.table["b"].to("g/cm2")
+        param_c = self.table["c"].to("km")
+
+        # build list of column density functions and their derivatives:
+        self._funcs = [
+            partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
+            for i, f in enumerate([_exponential] * 4 + [_linear])
+        ]
+        self._d_funcs = [
+            partial(f, a=param_a[i], b=param_b[i], c=param_c[i])
+            for i, f in enumerate([_d_exponential] * 4 + [_d_linear])
+        ]
+
+    @classmethod
+    def from_array(cls, array: np.ndarray, metadata: Dict = None):
+        """construct from a 5x5 array as provided by eventio"""
+
+        if metadata is None:
+            metadata = {}
+
+        if array.shape != (5, 5):
+            raise ValueError("expected ndarray with shape (5,5)")
+
+        table = Table(
+            array,
+            names=["height", "a", "b", "c", "1/c"],
+            units=["cm", "g/cm2", "g/cm2", "cm", "cm-1"],
+            meta=metadata,
+        )
+
+        table.meta.update(
+            dict(
+                TAB_VER=1,
+                TAB_TYPE="ctapipe.atmosphere.FiveLayerAtmosphereDensityProfile",
+            )
+        )
+        return cls(table)
+
+    @u.quantity_input(height=u.m)
+    def __call__(self, height) -> u.Quantity:
+        which_func = np.digitize(height, self.table["height"]) - 1
+        condlist = [which_func == i for i in range(5)]
+        return u.Quantity(
+            -1
+            * np.piecewise(
+                height,
+                condlist=condlist,
+                funclist=self._d_funcs,
+            )
+        ).to(u.g / u.cm**3)
+
+    @u.quantity_input(height=u.m)
+    def integral(self, height) -> u.Quantity:
+        which_func = np.digitize(height, self.table["height"]) - 1
+        condlist = [which_func == i for i in range(5)]
+        return u.Quantity(
+            np.piecewise(
+                x=height,
+                condlist=condlist,
+                funclist=self._funcs,
+            )
+        ).to(u.g / u.cm**2)
+
+    def __repr__(self):
+        return (
+            f"{self.__class__.__name__}(meta={self.table.meta}, rows={len(self.table)})"
+        )
diff --git a/ctapipe/instrument/atmosphere.py b/ctapipe/instrument/atmosphere.py
index e0f85c60..76a9e467 100644
--- a/ctapipe/instrument/atmosphere.py
+++ b/ctapipe/instrument/atmosphere.py
@@ -3,6 +3,7 @@ Functions to retrieve and interpolate atmosphere profiles.
 """
 import numpy as np
 from astropy.units import Quantity
+from astropy.utils.decorators import deprecated
 from scipy.interpolate import interp1d
 
 from ctapipe.utils import get_table_dataset
@@ -10,6 +11,9 @@ from ctapipe.utils import get_table_dataset
 __all__ = ["get_atmosphere_profile_table", "get_atmosphere_profile_functions"]
 
 
+@deprecated(
+    since="v0.16.0", message="use ctapipe.atmosphere profiles from EventSource instead"
+)
 def get_atmosphere_profile_table(atmosphere_name="paranal"):
     """
     Get an atmosphere profile table
@@ -30,6 +34,9 @@ def get_atmosphere_profile_table(atmosphere_name="paranal"):
     return table
 
 
+@deprecated(
+    since="v0.16.0", message="use ctapipe.atmosphere profiles from EventSource instead"
+)
 def get_atmosphere_profile_functions(atmosphere_name="paranal", with_units=True):
     """
     Gives atmospheric profile as a continuous function thickness(
diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py
index ecc295b8..14f33a22 100644
--- a/ctapipe/io/datawriter.py
+++ b/ctapipe/io/datawriter.py
@@ -23,6 +23,7 @@ from ..core.traits import Bool, CaselessStrEnum, Float, Int, Path, Unicode
 from ..instrument import SubarrayDescription
 from . import EventSource, HDF5TableWriter, TableWriter
 from . import metadata as meta
+from .astropy_helpers import write_table
 from .datalevels import DataLevel
 from .simteleventsource import SimTelEventSource
 from .tableio import FixedPointColumnTransform, TelListToMaskTransform
@@ -273,6 +274,9 @@ class DataWriter(Component):
         self._write_scheduling_and_observation_blocks()
         if self._is_simulation:
             self._write_simulation_configuration()
+            self._write_atmosphere_profile(
+                "/simulation/service/atmosphere_density_profile"
+            )
 
     def __enter__(self):
         return self
@@ -797,3 +801,31 @@ class DataWriter(Component):
             context_dict[key] = value
 
         meta.write_to_hdf5(context_dict, self._writer.h5file)
+
+    def _write_atmosphere_profile(self, path):
+        """
+        write atmosphere profiles if they are in a tabular format
+
+        Parameters
+        ----------
+        path: str
+            path in the HDF5 file where to place the profile
+
+        """
+
+        profile = self.event_source.atmosphere_density_profile
+
+        if profile:
+            if hasattr(profile, "table"):
+                write_table(
+                    table=profile.table,
+                    h5file=self._writer.h5file,
+                    path=path,
+                    append=False,
+                )
+            else:
+                self.logger.warning(
+                    f"The AtmosphereDensityProfile type '{profile.__class__.__name__}' "
+                    "is not serializable. No atmosphere profile will be stored in the "
+                    "output file"
+                )
diff --git a/ctapipe/io/eventsource.py b/ctapipe/io/eventsource.py
index 2ce3bf9f..1debd1ca 100644
--- a/ctapipe/io/eventsource.py
+++ b/ctapipe/io/eventsource.py
@@ -7,6 +7,8 @@ from typing import Dict, Generator, List, Tuple
 
 from traitlets.config.loader import LazyConfigValue
 
+from ctapipe.atmosphere import AtmosphereDensityProfile
+
 from ..containers import (
     ArrayEventContainer,
     ObservationBlockContainer,
@@ -280,6 +282,20 @@ class EventSource(Component):
         """
         return list(self.observation_blocks.keys())
 
+    @property
+    def atmosphere_density_profile(self) -> AtmosphereDensityProfile:
+        """atmosphere density profile that can be integrated to
+        convert between h_max and X_max.  This should correspond
+        either to what was used in a simualtion, or a measurment
+        for use with observed data.
+
+        Returns
+        -------
+        AtmosphereDensityProfile:
+           profile to be used
+        """
+        return None
+
     @abstractmethod
     def _generator(self) -> Generator[ArrayEventContainer, None, None]:
         """
diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py
index d92f7c9c..ccae1d76 100644
--- a/ctapipe/io/hdf5eventsource.py
+++ b/ctapipe/io/hdf5eventsource.py
@@ -8,6 +8,7 @@ import numpy as np
 import tables
 from astropy.utils.decorators import lazyproperty
 
+from ctapipe.atmosphere import AtmosphereDensityProfile
 from ctapipe.instrument.optics import FocalLengthKind
 
 from ..containers import (
@@ -42,6 +43,7 @@ from ..core import Container, Field
 from ..core.traits import UseEnum
 from ..instrument import SubarrayDescription
 from ..utils import IndexFinder
+from .astropy_helpers import read_table
 from .datalevels import DataLevel
 from .eventsource import EventSource
 from .hdf5tableio import HDF5TableReader, get_column_attrs
@@ -86,6 +88,33 @@ def get_hdf5_datalevels(h5file):
     return tuple(datalevels)
 
 
+def read_atmosphere_density_profile(
+    h5file: tables.File, path="/simulation/service/atmosphere_density_profile"
+):
+    """return a subclass of AtmosphereDensityProfile by
+    reading a table in a h5 file
+
+    Parameters
+    ----------
+    h5file: tables.File
+        file handle of HDF5 file
+    path: str
+        path in the file where the serialized model is stored as an
+        astropy.table.Table
+
+    Returns
+    -------
+    AtmosphereDensityProfile:
+        subclass depending on type of table
+    """
+
+    if path not in h5file:
+        return None
+
+    table = read_table(h5file=h5file, path=path)
+    return AtmosphereDensityProfile.from_table(table)
+
+
 class HDF5EventSource(EventSource):
     """
     Event source for files in the ctapipe DL1 format.
@@ -248,6 +277,10 @@ class HDF5EventSource(EventSource):
     def datalevels(self):
         return get_hdf5_datalevels(self.file_)
 
+    @lazyproperty
+    def atmosphere_density_profile(self) -> AtmosphereDensityProfile:
+        return read_atmosphere_density_profile(self.file_)
+
     @lazyproperty
     def obs_ids(self):
         return list(np.unique(self.file_.root.dl1.event.subarray.trigger.col("obs_id")))
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 4dd57cc2..ad7b35c5 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -1,17 +1,25 @@
 import enum
 import warnings
+from contextlib import nullcontext
+from enum import Enum, auto, unique
 from gzip import GzipFile
 from io import BufferedReader
 from pathlib import Path
-from typing import Dict
+from typing import Dict, Optional, Union
 
 import numpy as np
 from astropy import units as u
 from astropy.coordinates import Angle, EarthLocation
+from astropy.table import Table
 from astropy.time import Time
 from eventio.file_types import is_eventio
 from eventio.simtel.simtelfile import SimTelFile
 
+from ..atmosphere import (
+    AtmosphereDensityProfile,
+    FiveLayerAtmosphereDensityProfile,
+    TableAtmosphereDensityProfile,
+)
 from ..calib.camera.gainselection import GainSelector
 from ..containers import (
     ArrayEventContainer,
@@ -39,6 +47,7 @@ from ..containers import (
 )
 from ..coordinates import CameraFrame
 from ..core import Map
+from ..core.provenance import Provenance
 from ..core.traits import Bool, Float, Undefined, UseEnum, create_class_enum_trait
 from ..instrument import (
     CameraDescription,
@@ -61,7 +70,12 @@ from ..reco.impact_distance import shower_impact_distance
 from .datalevels import DataLevel
 from .eventsource import EventSource
 
-__all__ = ["SimTelEventSource"]
+__all__ = [
+    "SimTelEventSource",
+    "read_atmosphere_profile_from_simtel",
+    "AtmosphereProfileKind",
+    "MirrorClass",
+]
 
 # Mapping of SimTelArray Calibration trigger types to EventType:
 # from simtelarray: type Dark (0), pedestal (1), in-lid LED (2) or laser/LED (3+) data.
@@ -274,6 +288,113 @@ def apply_simtel_r1_calibration(
     return r1_waveforms, selected_gain_channel
 
 
+@unique
+class AtmosphereProfileKind(Enum):
+    """
+    choice for which kind of atmosphere density profile to load if more than
+    one is present"""
+
+    #: Don't load a profile
+    NONE = auto()
+
+    #: Use a TableAtmosphereDensityProfile
+    TABLE = auto()
+
+    #: Use a FiveLayerAtmosphereDensityProfile
+    FIVELAYER = auto()
+
+    #: Try TABLE first, and if it doesn't exist, use FIVELAYER
+    AUTO = auto()
+
+
+def read_atmosphere_profile_from_simtel(
+    simtelfile: Union[str, Path, SimTelFile], kind=AtmosphereProfileKind.AUTO
+) -> Optional[TableAtmosphereDensityProfile]:
+    """Read an atmosphere profile from a SimTelArray file as an astropy Table
+
+    Parameters
+    ----------
+    simtelfile: str | SimTelFile
+        filename of a SimTelArray file containing an atmosphere profile
+    kind: AtmosphereProfileKind | str
+        which type of model to load. In AUTO mode, table is tried first,
+        and if it doesn't exist, fivelayer is used.
+
+    Returns
+    -------
+    Optional[TableAtmosphereDensityProfile]:
+        Profile read from a table, with interpolation
+
+    """
+    if not isinstance(kind, AtmosphereProfileKind):
+        raise TypeError(
+            "read_atmosphere_profile_from_simtel: kind should be a AtmosphereProfileKind"
+        )
+
+    profiles = []
+
+    if kind == AtmosphereProfileKind.NONE:
+        return None
+
+    if isinstance(simtelfile, (str, Path)):
+        context_manager = SimTelFile(simtelfile)
+        Provenance().add_input_file(
+            filename=simtelfile, role="ctapipe.atmosphere.AtmosphereDensityProfile"
+        )
+
+    else:
+        context_manager = nullcontext(simtelfile)
+
+    with context_manager as simtel:
+
+        if (
+            not hasattr(simtel, "atmospheric_profiles")
+            or len(simtel.atmospheric_profiles) == 0
+        ):
+            return []
+
+        for atmo in simtel.atmospheric_profiles:
+
+            metadata = dict(
+                observation_level=atmo["obslevel"] * u.cm,
+                atmosphere_id=atmo["id"],
+                atmosphere_name=atmo["name"].decode("utf-8"),
+                atmosphere_height=atmo["htoa"] * u.cm,
+            )
+
+            if "altitude_km" in atmo and kind in {
+                AtmosphereProfileKind.TABLE,
+                AtmosphereProfileKind.AUTO,
+            }:
+                table = Table(
+                    dict(
+                        height=atmo["altitude_km"] * u.km,
+                        density=atmo["rho"] * u.g / u.cm**3,
+                        column_density=atmo["thickness"] * u.g / u.cm**2,
+                    ),
+                    meta=metadata,
+                )
+                profiles.append(TableAtmosphereDensityProfile(table=table))
+
+            elif (
+                kind == AtmosphereProfileKind.FIVELAYER
+                and "five_layer_atmosphere" in atmo
+            ):
+                profiles.append(
+                    FiveLayerAtmosphereDensityProfile.from_array(
+                        atmo["five_layer_atmosphere"], metadata=metadata
+                    )
+                )
+
+            else:
+                raise ValueError(f"Couldn't  load requested profile '{kind}'")
+
+    if profiles:
+        return profiles[0]
+
+    return None
+
+
 class SimTelEventSource(EventSource):
     """
     Read events from a SimTelArray data file (in EventIO format).
@@ -348,6 +469,17 @@ class SimTelEventSource(EventSource):
         ),
     ).tag(config=True)
 
+    atmosphere_profile_choice = UseEnum(
+        AtmosphereProfileKind,
+        default_value=AtmosphereProfileKind.AUTO,
+        help=(
+            "Which type of atmosphere density profile to load "
+            "from the file, in case more than one exists.  If set "
+            "to AUTO, TABLE will be attempted first and if missing, "
+            "FIVELAYER will be loaded."
+        ),
+    ).tag(config=True)
+
     def __init__(self, input_url=Undefined, config=None, parent=None, **kwargs):
         """
         EventSource for simtelarray files using the pyeventio library.
@@ -390,6 +522,11 @@ class SimTelEventSource(EventSource):
         self.gain_selector = GainSelector.from_name(
             self.gain_selector_type, parent=self
         )
+
+        self._atmosphere_density_profile = read_atmosphere_profile_from_simtel(
+            self.file_, kind=self.atmosphere_profile_choice
+        )
+
         self.log.debug(f"Using gain selector {self.gain_selector}")
 
     def __exit__(self, exc_type, exc_val, exc_tb):
@@ -410,6 +547,10 @@ class SimTelEventSource(EventSource):
     def simulation_config(self) -> Dict[int, SimulationConfigContainer]:
         return self._simulation_config
 
+    @property
+    def atmosphere_density_profile(self) -> AtmosphereDensityProfile:
+        return self._atmosphere_density_profile
+
     @property
     def observation_blocks(self) -> Dict[int, ObservationBlockContainer]:
         """
diff --git a/docs/conf.py b/docs/conf.py
index 2d507da0..bf22faa3 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -116,6 +116,7 @@ nitpick_ignore = [
     ("py:obj", "name"),
     ("py:class", "astropy.coordinates.baseframe.BaseCoordinateFrame"),
     ("py:class", "astropy.table.table.Table"),
+    ("py:class", "eventio.simtel.simtelfile.SimTelFile"),
 ]
 
 # The suffix(es) of source filenames.
diff --git a/docs/ctapipe_api/atmosphere/index.rst b/docs/ctapipe_api/atmosphere/index.rst
new file mode 100644
index 00000000..78ca8259
--- /dev/null
+++ b/docs/ctapipe_api/atmosphere/index.rst
@@ -0,0 +1,16 @@
+.. _atmosphere:
+
+==================================================
+Atmosphere Models (`~ctapipe.atmosphere`)
+==================================================
+
+.. currentmodule:: ctapipe.atmosphere
+
+Models of the atmosphere useful for transforming between *column densities*
+(atmosphere depths :math:`X`, in units of mass per area) and *heights* (in distance
+above sea-level units).
+
+Reference/API
+=============
+
+.. automodapi:: ctapipe.atmosphere


Test patch:
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index 80a23be0..04cf133d 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -1,3 +1,5 @@
+""" tests of SimTelEventSource """
+# pylint: disable=import-outside-toplevel
 import copy
 from itertools import zip_longest
 from pathlib import Path
@@ -13,7 +15,12 @@ from ctapipe.calib.camera.gainselection import ThresholdGainSelector
 from ctapipe.instrument.camera.geometry import UnknownPixelShapeWarning
 from ctapipe.instrument.optics import ReflectorShape
 from ctapipe.io import DataLevel
-from ctapipe.io.simteleventsource import SimTelEventSource, apply_simtel_r1_calibration
+from ctapipe.io.simteleventsource import (
+    AtmosphereProfileKind,
+    SimTelEventSource,
+    apply_simtel_r1_calibration,
+    read_atmosphere_profile_from_simtel,
+)
 from ctapipe.utils import get_dataset_path
 
 gamma_test_large_path = get_dataset_path("gamma_test_large.simtel.gz")
@@ -517,6 +524,45 @@ def test_simtel_no_metadata(monkeypatch):
         assert all([t.optics.name.startswith("UNKNOWN") for t in subarray.tel.values()])
 
 
+def test_load_atmosphere_from_simtel(prod5_gamma_simtel_path):
+    """
+    Load atmosphere from a SimTelEventSource
+    """
+    from ctapipe.atmosphere import (
+        FiveLayerAtmosphereDensityProfile,
+        TableAtmosphereDensityProfile,
+    )
+
+    profile = read_atmosphere_profile_from_simtel(
+        prod5_gamma_simtel_path, kind=AtmosphereProfileKind.AUTO
+    )
+    assert isinstance(profile, TableAtmosphereDensityProfile)
+
+    profile = read_atmosphere_profile_from_simtel(
+        prod5_gamma_simtel_path, kind=AtmosphereProfileKind.TABLE
+    )
+    assert isinstance(profile, TableAtmosphereDensityProfile)
+
+    profile = read_atmosphere_profile_from_simtel(
+        prod5_gamma_simtel_path, kind=AtmosphereProfileKind.FIVELAYER
+    )
+    assert isinstance(profile, FiveLayerAtmosphereDensityProfile)
+
+    # old simtel files don't have the profile in them, so a null list should be
+    # returned
+    simtel_path_old = get_dataset_path("gamma_test_large.simtel.gz")
+    profile = read_atmosphere_profile_from_simtel(simtel_path_old)
+    assert not profile
+
+
+def test_atmosphere_profile(prod5_gamma_simtel_path):
+    """check that for a file with a profile in it that we get it back"""
+    from ctapipe.atmosphere import AtmosphereDensityProfile
+
+    with SimTelEventSource(prod5_gamma_simtel_path) as source:
+        assert isinstance(source.atmosphere_density_profile, AtmosphereDensityProfile)
+
+
 @pytest.mark.parametrize("sign", (-1, 1))
 def test_float32_pihalf(sign):
     float32_pihalf = np.float32(sign * np.pi / 2)
diff --git a/ctapipe/tests/test_atmosphere.py b/ctapipe/tests/test_atmosphere.py
new file mode 100644
index 00000000..2fc06484
--- /dev/null
+++ b/ctapipe/tests/test_atmosphere.py
@@ -0,0 +1,135 @@
+"""
+Checks for the atmosphere module
+"""
+
+# pylint: disable=import-outside-toplevel
+from unittest.mock import patch
+
+import astropy.units as u
+import numpy as np
+import pytest
+
+import ctapipe.atmosphere as atmo
+from ctapipe.utils import get_dataset_path
+
+SIMTEL_PATH = get_dataset_path(
+    "gamma_20deg_0deg_run2___cta-prod5-paranal_desert"
+    "-2147m-Paranal-dark_cone10-100evts.simtel.zst"
+)
+
+
+def get_simtel_profile_from_eventsource():
+    """get a TableAtmosphereDensityProfile from a simtel file"""
+    from ctapipe.io import EventSource
+
+    with EventSource(SIMTEL_PATH) as source:
+        return source.atmosphere_density_profile
+
+
+@pytest.fixture(scope="session", name="table_profile")
+def fixture_table_profile():
+    """a table profile for testing"""
+    return get_simtel_profile_from_eventsource()
+
+
+def get_simtel_fivelayer_profile():
+    """
+    get a sample 3-layer profile
+    """
+    from ctapipe.io.simteleventsource import (
+        AtmosphereProfileKind,
+        read_atmosphere_profile_from_simtel,
+    )
+
+    return read_atmosphere_profile_from_simtel(
+        SIMTEL_PATH, kind=AtmosphereProfileKind.FIVELAYER
+    )
+
+
+@pytest.mark.parametrize(
+    "density_model",
+    [
+        atmo.ExponentialAtmosphereDensityProfile(),
+        get_simtel_profile_from_eventsource(),
+        get_simtel_fivelayer_profile(),
+    ],
+)
+def test_models(density_model):
+    """check that a set of model classes work"""
+
+    # test we can convert to correct units
+    density_model(10 * u.km).to(u.kg / u.m**3)
+
+    # ensure units are properly taken into account
+    assert np.isclose(density_model(1 * u.km), density_model(1000 * u.m))
+
+    # check we can also compute the integral
+    column_density = density_model.integral(10 * u.km)
+    assert column_density.unit.is_equivalent(u.g / u.cm**2)
+
+    assert np.isclose(
+        density_model.integral(1 * u.km), density_model.integral(1000 * u.m)
+    )
+
+    with patch("matplotlib.pyplot.show"):
+        density_model.peek()
+
+
+def test_exponential_model():
+    """check exponential models"""
+
+    density_model = atmo.ExponentialAtmosphereDensityProfile(
+        scale_height=10 * u.m, scale_density=0.00125 * u.g / u.cm**3
+    )
+    assert np.isclose(density_model(1_000_000 * u.km), 0 * u.g / u.cm**3)
+    assert np.isclose(density_model(0 * u.km), density_model.scale_density)
+
+
+def test_table_model_interpolation(table_profile):
+    """check that interpolation is reasonable"""
+
+    np.testing.assert_allclose(
+        table_profile(table_profile.table["height"].to("km")),
+        table_profile.table["density"].to("g cm-3"),
+    )
+
+    # check that fine interpolation up to 100 km :
+    height_fine = np.linspace(0, 100, 1000) * u.km
+    assert np.isfinite(table_profile.integral(height_fine)).all()
+
+
+def test_against_reference():
+    """
+    Test five-layer and table methods against a reference analysis from
+    SimTelArray.  Data from communication with K. Bernloehr.
+
+    See https://github.com/cta-observatory/ctapipe/pull/2000
+    """
+    from ctapipe.utils import get_table_dataset
+
+    reference_table = get_table_dataset(
+        "atmosphere_profile_comparison_from_simtelarray"
+    )
+
+    fit_reference = np.array(
+        [
+            [0.00 * 100000, -140.508, 1178.05, 994186, 0],
+            [9.75 * 100000, -18.4377, 1265.08, 708915, 0],
+            [19.00 * 100000, 0.217565, 1349.22, 636143, 0],
+            [46.00 * 100000, -0.000201796, 703.745, 721128, 0],
+            [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],
+        ]
+    )
+
+    profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)
+
+    height = reference_table["Altitude_km"].to("km")
+
+    np.testing.assert_allclose(
+        1.0 - profile_5(height) / reference_table["rho_5"], 0, atol=1e-5
+    )
+    np.testing.assert_allclose(
+        1.0 - profile_5.line_of_sight_integral(height) / reference_table["thick_5"],
+        0,
+        atol=1e-5,
+    )