Repository: kyleaoman/martini (language: python)

Issue:
Support FIRE
FIRE has a public data release and there seems to be some interest in using MARTINI, time to add a new source module...


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 0d46d29..fa35aa2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,6 +9,7 @@
 /docs/build/
 .ipynb_checkpoints
 examples/RefL0025N0376/*
+examples/m11e_res7100/*
 examples/*.fits
 examples/*.FITS
 examples/*.hdf5
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 40a8a6f..3a737f5 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -9,7 +9,7 @@ version: 2
 build:
   os: ubuntu-22.04
   tools:
-    python: "3.10"
+    python: "3.12"
 
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
diff --git a/README.rst b/README.rst
index 7e94719..f0fab01 100644
--- a/README.rst
+++ b/README.rst
@@ -136,7 +136,7 @@ Installation Notes
 
 .. INSTALLATION_NOTES_START_LABEL
 
-MARTINI works with ``python3`` (version ``3.8`` or higher), and does not support ``python2``.
+MARTINI works with ``python3`` (version ``3.9`` or higher), and does not support ``python2``.
 
 Stable releases are available via PyPI_:
 
@@ -149,7 +149,7 @@ and the numbered releases (starting from ``2.0.0``) on github. The github main b
 .. _PyPI: https://pypi.org/project/astromartini/
 .. _kyle.a.oman@durham.ac.uk: mailto:kyle.a.oman@durham.ac.uk
 
-The easiest way to install MARTINI is from PyPI by doing ``python3 -m pip install astromartini``. Output to ``.fits`` files is supported by default; if output to ``.hdf5`` format is desired use ``python3 -m pip install "astromartini[hdf5_output]"`` instead. This will also handle the installation of the required dependencies. Other optional features require additional dependencies hosted on PyPI. In particular, EAGLE, Illustris/TNG and Simba users who wish to use the custom source modules for those simulations in MARTINI can automatically install the optional dependencies with ``python3 -m pip install "astromartini[eaglesource]"``, or ``python3 -m pip install "astromartini[simbasource]"`` or ``python3 -m pip install "astromartini[tngsource]"``.
+The easiest way to install MARTINI is from PyPI by doing ``python3 -m pip install astromartini``. Output to ``.fits`` files is supported by default; if output to ``.hdf5`` format is desired use ``python3 -m pip install "astromartini[hdf5_output]"`` instead. This will also handle the installation of the required dependencies. Other optional features require additional dependencies hosted on PyPI. In particular, EAGLE, Illustris/TNG, Simba and FIRE users who wish to use the custom source modules for those simulations in MARTINI can automatically install the optional dependencies with ``python3 -m pip install "astromartini[eaglesource]"``, ``python3 -m pip install "astromartini[simbasource]"``, ``python3 -m pip install "astromartini[tngsource]"``, or ``python3 -m pip install "astromartini[firesource]"``.
 
 .. INSTALLATION_NOTES_END_LABEL
 
@@ -192,15 +192,13 @@ See the help for |martini.Martini| for an example script to configure MARTINI an
 
 MARTINI has (so far) been successfully run on the output of these simulations:
 
-* EAGLE
-* APOSTLE
-* C-EAGLE/Hydrangea
-* Illustris
-* IllustrisTNG
-* Auriga
-* MaGICC (and therefore in principle NIHAO)
-* Magneticum
+* EAGLE (also APOSTLE, C-EAGLE/Hydrangea)
+* IllustrisTNG (also Illustris, Auriga)
 * Simba
+* FIRE
+* Magneticum
+* MaGICC (and therefore in principle NIHAO)
+* Colibre (test runs)
 
 I attempt to support publicly available simulations with a customized source module. If your simulation is public and not supported, please `contact me`_. Currently custom source modules exist for:
 
@@ -208,10 +206,12 @@ I attempt to support publicly available simulations with a customized source mod
 
 * EAGLE (|martini.sources.EAGLESource|)
 * IllustrisTNG (|martini.sources.TNGSource|; also works with Illustris)
-* Magneticum (|martini.sources.MagneticumSource|)
+* FIRE (|martini.sources.FIRESource|)
 * Simba (|martini.sources.SimbaSource|)
+* Magneticum (|martini.sources.MagneticumSource|)
+* Colibre (|martini.sources.ColibreSource|)
 
-Example notebooks_ are available for supported simulations.
+Example notebooks_ are available for supported, publicly available simulations.
 
 .. _notebooks: https://github.com/kyleaoman/martini/tree/main/examples
 
@@ -220,6 +220,8 @@ Example notebooks_ are available for supported simulations.
 .. |martini.Martini| replace:: `martini.Martini <https://martini.readthedocs.io/en/latest/modules/martini.martini.html#martini.martini.Martini>`__
 .. |martini.sources.EAGLESource| replace:: `martini.sources.EAGLESource <https://martini.readthedocs.io/en/latest/modules/martini.sources.eagle_source.html#martini.sources.eagle_source.EAGLESource>`__
 .. |martini.sources.TNGSource| replace:: `martini.sources.TNGSource <https://martini.readthedocs.io/en/latest/modules/martini.sources.tng_source.html#martini.sources.tng_source.TNGSource>`__
+.. |martini.sources.FIRESource| replace:: `martini.sources.FIRESource <https://martini.readthedocs.io/en/latest/modules/martini.sources.fire_source.html#martini.sources.fire_source.FIRESource>`__					 
 .. |martini.sources.MagneticumSource| replace:: `martini.sources.MagneticumSource <https://martini.readthedocs.io/en/latest/modules/martini.sources.magneticum_source.html#martini.sources.magneticum_source.MagneticumSource>`__
 .. |martini.sources.SimbaSource| replace:: `martini.sources.SimbaSource <https://martini.readthedocs.io/en/latest/modules/martini.sources.simba_source.html#martini.sources.simba_source.SimbaSource>`__
 .. |martini.sources.SOSource| replace:: `martini.sources.SOSource <https://martini.readthedocs.io/en/latest/modules/martini.sources.so_source.html#martini.sources.so_source.SOSource>`__
+.. |martini.sources.ColibreSource| replace:: `martini.sources.ColibreSource <https://martini.readthedocs.io/en/latest/modules/martini.sources.colibre_source.html#martini.sources.colibre_source.ColibreSource>`__
diff --git a/codecov.yml b/codecov.yml
index 022300f..d3b094c 100644
--- a/codecov.yml
+++ b/codecov.yml
@@ -10,4 +10,5 @@ ignore:
   - "martini/sources/_illustris_tools.py"
   - "martini/sources/tng_source.py"
   - "martini/sources/eagle_source.py"
-  - "martini/sources/simba_source.py"
\ No newline at end of file
+  - "martini/sources/simba_source.py"
+  - "martini/sources/fire_source.py"
diff --git a/docs/source/getting_started/index.rst b/docs/source/getting_started/index.rst
index 5b76783..9243f29 100644
--- a/docs/source/getting_started/index.rst
+++ b/docs/source/getting_started/index.rst
@@ -29,3 +29,5 @@ Quickstart guide
 .. |martini.sources.MagneticumSource| replace:: :class:`~martini.sources.magneticum_source.MagneticumSource`
 .. |martini.sources.SimbaSource| replace:: :class:`~martini.sources.simba_source.SimbaSource`
 .. |martini.sources.SOSource| replace:: :class:`~martini.sources.so_source.SOSource`
+.. |martini.sources.FIRESource| replace:: :class:`~martini.sources.fire_source.FIRESource`
+.. |martini.sources.ColibreSource| replace:: :class:`~martini.sources.colibre_source.ColibreSource`
diff --git a/docs/source/modules/martini.sources.fire_source.rst b/docs/source/modules/martini.sources.fire_source.rst
new file mode 100644
index 0000000..7411e54
--- /dev/null
+++ b/docs/source/modules/martini.sources.fire_source.rst
@@ -0,0 +1,7 @@
+fire_source module
+===================
+
+.. automodule:: martini.sources.fire_source
+  :members:
+  :show-inheritance:
+  :inherited-members:
diff --git a/docs/source/modules/martini.sources.rst b/docs/source/modules/martini.sources.rst
index 0efe50b..c120346 100644
--- a/docs/source/modules/martini.sources.rst
+++ b/docs/source/modules/martini.sources.rst
@@ -9,6 +9,7 @@ sources module
    martini.sources.swiftgalaxy_source
    martini.sources.colibre_source
    martini.sources.eagle_source
+   martini.soures.fire_source
    martini.sources.magneticum_source
    martini.sources.simba_source
    martini.sources.tng_source
diff --git a/docs/source/sources/index.rst b/docs/source/sources/index.rst
index 9534745..7ee90b3 100644
--- a/docs/source/sources/index.rst
+++ b/docs/source/sources/index.rst
@@ -6,11 +6,12 @@ In the context of mock observations of a simulation, the source is the collectio
 Sources in MARTINI
 ------------------
 
-MARTINI has a collection of classes used to store and manipulate the properties of the source particles. All of these inherit from the :class:`martini.sources.sph_source.SPHSource` class. This base class can be created by providing arrays containing the particle coordinates, velocities, masses, etc. The other source classes are tailored to specific simulations or data formats and will typically expect one or more filenames as input, along with an identifier specifying a group of particles (as defined in a group catalogue, for example). These classes will then take care of reading data from the input files, annotating them with units, doing any required calculations or conversions, and so on. There are source classes for the EAGLE_, IllustrisTNG_, Simba_, Colibre and Magneticum_ simulations. The SWIFTGalaxy_ interface to `SWIFT`_-based simulations with a variety of halo catalogue formats is also supported.
+MARTINI has a collection of classes used to store and manipulate the properties of the source particles. All of these inherit from the :class:`martini.sources.sph_source.SPHSource` class. This base class can be created by providing arrays containing the particle coordinates, velocities, masses, etc. The other source classes are tailored to specific simulations or data formats and will typically expect one or more filenames as input, along with an identifier specifying a group of particles (as defined in a group catalogue, for example). These classes will then take care of reading data from the input files, annotating them with units, doing any required calculations or conversions, and so on. There are source classes for the EAGLE_, IllustrisTNG_, Simba_, FIRE_, Colibre and Magneticum_ simulations. The SWIFTGalaxy_ interface to `SWIFT`_-based simulations with a variety of halo catalogue formats is also supported.
 
 .. _EAGLE: https://icc.dur.ac.uk/Eagle
 .. _IllustrisTNG: https://www.tng-project.org
 .. _Simba: http://simba.roe.ac.uk
+.. _FIRE: https://fire.northwestern.edu/
 .. _Magneticum: http://www.magneticum.org/
 .. _SWIFTGalaxy: https://github.com/SWIFTSIM/swiftgalaxy
 .. _SWIFT: https://github.com/SWIFTSIM/SWIFT
@@ -69,11 +70,11 @@ The smoothing lengths are defined as the full width at half-maximum (FWHM) of th
 Simulation-specific source modules
 ++++++++++++++++++++++++++++++++++
 
-MARTINI provides source modules to simplify working with publicly-available simulation data sets. These currently include the :class:`~martini.sources.eagle_source.EAGLESource`, :class:`~martini.sources.tng_source.TNGSource` and :class:`~martini.sources.simba_source.SimbaSource`. Example usage including how to obtain publicly available data is provided in a set of `Jupyter notebooks`_.
+MARTINI provides source modules to simplify working with publicly-available simulation data sets. These currently include the :class:`~martini.sources.eagle_source.EAGLESource`, :class:`~martini.sources.tng_source.TNGSource`, :class:`~martini.sources.simba_source.SimbaSource` and :class:`~martini.sources.fire_source.FIRESource`. Example usage including how to obtain publicly available data is provided in a set of `Jupyter notebooks`_.
 
 The Magneticum_ simulations are supported with the :class:`~martini.sources.magneticum_source.MagneticumSource`, however detailed examples are not provided due to the lack of publicly available snapshot data. Interested users can refer to the API documentation of the class.
 
-There is also a :class:`~martini.sources.swiftgalaxy_source.SWIFTGalaxySource` source module for SWIFT_ simulation data in conjunction with Velociraptor_ or Caesar_ halo catalogues; support for `HBT+`_ and SOAP_ is also planned. A :class:`~martini.sources.colibre_source.ColibreSource` module exists, but since this galaxy formation model is still under development this is currently experimental functionality.
+There is also a :class:`~martini.sources.swiftgalaxy_source.SWIFTGalaxySource` source module for SWIFT_ simulation data in conjunction with Velociraptor_ or Caesar_ halo catalogues; support for `HBT+`_ and/or SOAP_ is also planned. A :class:`~martini.sources.colibre_source.ColibreSource` module exists, but since this galaxy formation model is still under development this is currently experimental functionality.
 
 Finally, there is an :class:`~martini.sources.so_source.SOSource` module that interfaces with the SimObj_ package. That package is no longer maintained, but may facilitate working with some simulations including APOSTLE, C-EAGLE/Hydrangea and Auriga.
 
diff --git a/examples/README.rst b/examples/README.rst
index 674d3af..3b7d8e8 100644
--- a/examples/README.rst
+++ b/examples/README.rst
@@ -8,7 +8,9 @@ You can view the example notebooks online using nbviewer_.
 - EAGLE_
 - Simba_
 - TNG_
+- FIRE_
 
 .. _EAGLE: https://nbviewer.jupyter.org/github/kyleaoman/martini/blob/main/examples/martini_eagle.ipynb
 .. _Simba: https://nbviewer.jupyter.org/github/kyleaoman/martini/blob/main/examples/martini_simba.ipynb
 .. _TNG: https://nbviewer.jupyter.org/github/kyleaoman/martini/blob/main/examples/martini_TNG.ipynb
+.. _FIRE: https://nbviewer.jupyter.org/github/kyleaoman/martini/blob/main/examples/martini_fire. ipynb
diff --git a/examples/martini_TNG.ipynb b/examples/martini_TNG.ipynb
index 7dd93d0..0cbfce6 100644
--- a/examples/martini_TNG.ipynb
+++ b/examples/martini_TNG.ipynb
@@ -763,7 +763,7 @@
     "plt.colorbar(im1, ax=sp1, label=f\"mom0 [{mom0_unit}]\")\n",
     "im2 = sp2.imshow(\n",
     "    (mom1 * mask).to_value(mom1_unit),\n",
-    "    cmap=\"RdBu\",\n",
+    "    cmap=\"RdBu_r\",\n",
     "    vmin=-np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
     "    vmax=np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
     ")\n",
diff --git a/examples/martini_eagle.ipynb b/examples/martini_eagle.ipynb
index cc39d60..54e500a 100644
--- a/examples/martini_eagle.ipynb
+++ b/examples/martini_eagle.ipynb
@@ -863,7 +863,7 @@
     "plt.colorbar(im1, ax=sp1, label=f\"mom0 [{mom0_unit}]\")\n",
     "im2 = sp2.imshow(\n",
     "    (mom1 * mask).to_value(mom1_unit),\n",
-    "    cmap=\"RdBu\",\n",
+    "    cmap=\"RdBu_r\",\n",
     "    vmin=-np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
     "    vmax=np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
     ")\n",
diff --git a/examples/martini_fire.ipynb b/examples/martini_fire.ipynb
new file mode 100644
index 0000000..2eedc9e
--- /dev/null
+++ b/examples/martini_fire.ipynb
@@ -0,0 +1,1089 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this short tutorial we will install and use [MARTINI](https://martini.readthedocs.io/en/latest/), an analysis package for creating mock HI-data cubes similar to radio interferometer data, written by Kyle Oman (kyle.a.oman@durham.ac.uk). This example uses input from the [FIRE-2](https://ui.adsabs.harvard.edu/abs/2018MNRAS.480..800H/abstract) simulations. The data are publicly available and hosted at [flathub.flatironinstitute.org/fire](https://flathub.flatironinstitute.org/fire)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "![MARTINI](https://github.com/kyleaoman/martini/raw/main/martini_banner.png)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "MARTINI is a modular package for the creation of synthetic resolved HI line observations (data cubes) of smoothed-particle hydrodynamics simulations of galaxies. The various aspects of the mock-observing process are divided logically into sub-modules handling the data cube, source, beam, noise, spectral model and SPH kernel. MARTINI is object-oriented: each sub-module provides a class (or classes) which can be configured as desired. For most sub-modules, base classes are provided to allow for straightforward customization. Instances of each sub-module class are given as parameters to the Martini class; a mock observation is then constructed by calling a handful of functions to execute the desired steps in the mock-observing process.\n",
+    "\n",
+    "This tutorial focuses on particulars related to working with the FIRE simulations. More general information is available in the MARTINI documentation, [hosted on ReadTheDocs](https://martini.readthedocs.io/en/latest/)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Installation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "MARTINI requires `python3` version `3.8` or higher.\n",
+    "\n",
+    "The following command will use `pip` to download and install [MARTINI from pypi](https://pypi.org/project/astromartini/):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "skip-execution"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "\n",
+    "!{sys.executable} -m pip install \"astromartini[firesource]==2.1.3\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If you do not have superuser permissions or use a virtual environment, you may wish to add the --user flag.\n",
+    "With this command required dependencies will be fetched and installed automatically. Watch for error messages during installation. For greater control you may also install the dependencies by hand. These are: numpy, astropy, scipy and h5py.\n",
+    "\n",
+    "We'll also install a few other packages used in this notebook:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "skip-execution"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "!{sys.executable} -m pip install requests matplotlib beautifulsoup4"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This cell may be needed in some cases to display figures below:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's check that we can `import martini`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import martini"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If this produces errors, you may need to restart the Python kernel of this notebook so that it sees the recently installed packages (Kernel -> Restart in the menubar).\n",
+    "\n",
+    "We can run MARTINI's built-in demo to check that all of the basic functionality works:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "skip-execution"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from martini import demo\n",
+    "\n",
+    "demo()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "When run successfully, this will make a mock observation of a very simple analytic disc model and write some output to the working directory. Rather than inspect this toy model, let's look at a \"real\" simulation..."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## FIRE Data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This example uses data from the FIRE-2 simulations. These data are hosted at [https://flathub.flatironinstitute.org/fire](https://flathub.flatironinstitute.org/fire).\n",
+    "\n",
+    "In this example, we will use the [`m11e_res7100`](https://users.flatironinstitute.org/~chayward/fire2_public_release/core/m11e_res7100) simulation (80 GB). The following cell will download it directly, but note that a [Globus](https://www.globus.org/) endpoint is also provided (see [FlatHUB page](https://flathub.flatironinstitute.org/fire) for details) that can make transfers substantially faster:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "import urllib\n",
+    "import re\n",
+    "import requests\n",
+    "from tqdm import tqdm\n",
+    "import bs4\n",
+    "\n",
+    "data_url = \"https://users.flatironinstitute.org/~chayward/fire2_public_release/core/\"\n",
+    "simulation_directory = \"m11e_res7100/\"\n",
+    "\n",
+    "\n",
+    "def direct_download(url, fname):\n",
+    "    with requests.get(url) as r:\n",
+    "        r.raise_for_status()\n",
+    "        with open(fname, \"wb\") as f:\n",
+    "            f.write(r.content)\n",
+    "\n",
+    "\n",
+    "def chunked_download(url, fname):\n",
+    "    with requests.get(url, stream=True) as r:\n",
+    "        r.raise_for_status()\n",
+    "        total_size_in_bytes = int(r.headers.get(\"content-length\", 0))\n",
+    "        chunk_size = 8192\n",
+    "        progress_bar = tqdm(total=total_size_in_bytes, unit=\"iB\", unit_scale=True)\n",
+    "        with open(fname, \"wb\") as f:\n",
+    "            for chunk in r.iter_content(chunk_size=chunk_size):\n",
+    "                progress_bar.update(len(chunk))\n",
+    "                f.write(chunk)\n",
+    "    return\n",
+    "\n",
+    "\n",
+    "def directory_download(url_base, directory, chunked_if_larger_than=1024 * 1024 * 10):\n",
+    "    dir_url = urllib.parse.urljoin(url_base, directory)\n",
+    "    print(directory)\n",
+    "    os.mkdir(directory)\n",
+    "    r = requests.get(dir_url)\n",
+    "    data = bs4.BeautifulSoup(r.text, \"html.parser\")\n",
+    "    for link in data.find_all(\"a\"):\n",
+    "        link_url = urllib.parse.urljoin(dir_url, link[\"href\"])\n",
+    "        if \"..\" in link[\"href\"]:  # parent directory, skip\n",
+    "            continue\n",
+    "        elif \".\" in link[\"href\"]:  # a file, download\n",
+    "            target_file = urllib.parse.urljoin(directory, link[\"href\"])\n",
+    "            with requests.head(link_url) as h:\n",
+    "                h.raise_for_status()\n",
+    "                total_size_in_bytes = int(h.headers.get(\"content-length\", 0))\n",
+    "                print(target_file)\n",
+    "                if total_size_in_bytes > chunked_if_larger_than:\n",
+    "                    chunked_download(link_url, target_file)\n",
+    "                else:\n",
+    "                    direct_download(link_url, target_file)\n",
+    "        else:  # a directory, recurse\n",
+    "            directory_download(url_base, urllib.parse.urljoin(directory, link[\"href\"]))\n",
+    "\n",
+    "\n",
+    "if os.path.isdir(simulation_directory):\n",
+    "    print(\"Found directory matching simulation_directory locally, skipping download.\")\n",
+    "else:\n",
+    "    directory_download(data_url, simulation_directory)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If you have downloaded the data manually, edit this cell to specify the location of the base directory:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "simulation_directory = (\n",
+    "    simulation_directory  # e.g. simulation_directory = \"/path/to/m11e_res7100/\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The i/o tools for FIRE simulations allow a few ways of defining a snapshot, by its `\"index\"`, `\"redshift\"` or `\"scalefactor\"`. Here we'll use the `\"redshift\"` method and specify redshift `0`; the i/o tools will automatically select the closest snapshot to z=0. Since the FIRE suite is a set of zoom simulations (at least those in the FIRE-2 [public data release](https://flathub.flatironinstitute.org/fire)), the galaxy of interest is usually the first in the catalogue, corresponding to `host_number = 1`. Subsequent galaxies in the catalogue can be selected by increasing this number."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "snapshot = (\"redshift\", 0)\n",
+    "host_number = 1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## FIRE Example"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First, import some modules from Martini, and the units module from astropy."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from martini.sources import FIRESource\n",
+    "from martini import DataCube, Martini\n",
+    "from martini.beams import GaussianBeam\n",
+    "from martini.noise import GaussianNoise\n",
+    "from martini.spectral_models import GaussianSpectrum\n",
+    "from martini.sph_kernels import CubicSplineKernel\n",
+    "import astropy.units as U"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The different martini sub-modules need to be initialized, \n",
+    "see the [full documentation](https://martini.readthedocs.io/en/latest/) for details of all configuration options. A few suggested best-practices specific to FIRE are outlined below."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### SOURCE"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `FIRESource` module is built around the i/o tools for the FIRE simulations, which are built assuming that the target object in a zoom simulation is an object of interest. The default (and only) behaviour is to read in all of the particles of a given type in the snapshot. Initializing the source module will therefore have a memory usage similar to the size of the gas particle arrays in the snapshot (in this case about ?? GB). The memory footprint can be reduced by reducing precision to 32-bit floats with `convert_float32=True`. The FIRE i/o tools can print some diagnostic/progress messages, enable this with ``gizmo_io_verbose=True`` if desired.\n",
+    "\n",
+    "Galaxy positions and velocities are assigned by default using a \"shrinking spheres\" algorithm weighted by the particle masses (`assign_hosts=\"mass\"`), but other schemes are also available (potential-weighted with `assign_hosts=\"potential\"`, metallicity-weighted with `assign_hosts=\"massfraction.metals\"`, based on smoothed position across snapshots with `assign_hosts=\"track\"`, read directly from the Rockstar halo finder catalogues with `assign_hosts=\"halo\"`; setting `assign_hosts=True` will attempt to find a sensible default by first trying `\"track\"`, then `\"mass\"`).\n",
+    "\n",
+    "The `FIRESource` assumes the molecular gas partitioning scheme described in [Orr et al. (2018)](https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.3653O/abstract), Appendix B."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "source = FIRESource(\n",
+    "    simulation_directory=simulation_directory,\n",
+    "    snapshot=snapshot,\n",
+    "    host_number=host_number,\n",
+    "    assign_hosts=\"mass\",\n",
+    "    convert_float32=False,\n",
+    "    gizmo_io_verbose=False,\n",
+    "    distance=5.0 * U.Mpc,\n",
+    "    vpeculiar=0 * U.km / U.s,\n",
+    "    rotation={\"rotmat\": np.eye(3)},\n",
+    "    ra=0.0 * U.deg,\n",
+    "    dec=0.0 * U.deg,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The rotation argument above has been set to the identity matrix, so the source has the (random) orientation that it has within the simulation volume. The source class includes a function to make a quick plot to get an idea of the source's appearance:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "preview_fig_unrotated = source.preview(title=\"unrotated\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The preview function defaults to apertures in position and velocity that enclose all particles in the source (and therefore the entire zoom simulation), so this preview emphasizes the diffuse circumgalactic and intergalactic gas. The apertures can be set manually using the `lim` and `vlim` keywords to set the maximum absolute offsets in position and velocity relative to the source centre to be plotted. For example, restricting the aperture to 50kpc and 300km/s makes the disc more clearly visible."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "preview_fig_unrotated_zoom = source.preview(\n",
+    "    title=\"unrotated, zoomed-in\",\n",
+    "    lim=50 * U.kpc,\n",
+    "    vlim=300 * U.km / U.s,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This randomly-oriented viewing angle seems to be moderately inclined with respect to the disc. The source can be rotated to a different orientation. MARTINI's tool for quick/approximate manipulation of the orientation of the source aligns the source based on its angular momentum vector (\"L\"), for example:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "source.rotate(L_coords=(60 * U.deg, 90 * U.deg))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The rotation configuration takes an inclination (here 60deg) and rotation about the pole (here 90deg, relative to an arbitrary reference direction). The code attempts to\n",
+    "automatically align the galactic disk in the y-z plane by aligning\n",
+    "the angular momentum along the x-axis. The polar rotation is then\n",
+    "applied, and finally the disc inclined by a rotation around the\n",
+    "y-axis (the line of sight is along the x-axis). The automatic\n",
+    "alignment will work for typical reasonably isolated discs (no other galaxies with HI in the zoom region), but will\n",
+    "struggle when companions are present, when the angular momentum axis\n",
+    "is a poor tracer of the disc plane, and especially for satellites or galaxies in pairs."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "preview_fig_rotated_zoomed_in = source.preview(\n",
+    "    title=\"rotated, zoomed-in\",\n",
+    "    lim=50 * U.kpc,\n",
+    "    vlim=300 * U.km / U.s,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If finer control of the orientation is needed, derive the transformation from the simulation box coordinates (see [the documentation](https://martini.readthedocs.io/en/latest/) for examples) to the desired coordinates for the 'observation', keeping in mind that the line of sight is along the x-axis. This rotation matrix can then be passed to the rotate function as `rotmat=np.eye(3)` (here the identity rotation matrix used as a place holder). The rotation can also be provided when the source is initialized by using the `rotation` keyword argument.\n",
+    "\n",
+    "A common problem is deriving the inverse transform instead of the forward transform, if unexpected results are obtained, first try passing the transpose of the rotation matrix."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### DATACUBE"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It is usually advisable to set the centre of the cube to track the\n",
+    "centre of the source, as illustrated below. Note that the source\n",
+    "systemic velocity is set according to the distance, peculiar velocity, and Hubble's law.\n",
+    "These values can instead be set explicitly, if desired. A datacube\n",
+    "with 128x128 pixels usually takes a few minutes, depending on the number of particles. 1024x1024 can take\n",
+    "several hours. The number of channels has less influence on the\n",
+    "runtime. Most of the runtime is spent when `M.insert_source_in_cube` is\n",
+    "called below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "datacube = DataCube(\n",
+    "    n_px_x=384,\n",
+    "    n_px_y=384,\n",
+    "    n_channels=50,\n",
+    "    px_size=10.0 * U.arcsec,\n",
+    "    channel_width=5.0 * U.km * U.s**-1,\n",
+    "    velocity_centre=source.vsys,\n",
+    "    ra=source.ra,\n",
+    "    dec=source.dec,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### BEAM"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It is usually advisable to set the beam size to be ~3x the pixel\n",
+    "size. Note that the data cube is padded according to the size of the\n",
+    "beam, this usually results in the number of pixel rows printed in the\n",
+    "progress messages to differ from the requested dimensions. The\n",
+    "padding is required for accurate convolution with the beam, but\n",
+    "contains incorrect values after convolution and is discarded to\n",
+    "produce the final data cube of the requested size."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "beam = GaussianBeam(\n",
+    "    bmaj=30.0 * U.arcsec,\n",
+    "    bmin=30.0 * U.arcsec,\n",
+    "    bpa=0.0 * U.deg,\n",
+    "    truncate=3.0,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "### NOISE"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The noise is normally added before convolution with the beam (as\n",
+    "below in this example). The rms value passed is that corresponding to the desired noise level in the final data cube, in Jy/beam or equivalent units. To obtain consistent random realisations each time the code is run, we provide a random seed (integer)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "noise = GaussianNoise(\n",
+    "    rms=3.0e-8 * U.Jy * U.beam**-1,\n",
+    "    seed=0,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "### SPECTRAL MODEL"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `thermal` mode estimates the HI line width for each particle based on its properties (temperature, etc.). The 'subgrid' velocity dispersion can also be fixed to a constant value, e.g. `sigma=7 * U.km / U.s`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "spectral_model = GaussianSpectrum(sigma=\"thermal\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The calculation of the spectra (that will happen when the `Martini` module is initialized below) can be done in parallel by providing a keyword argument `ncpu=N`, where `N` is the number of CPUs to use. However, the details of the implementation mean that for small numbers of particles running in parallel tends to slow down the calculation, so turning this on should be done with care. Significant speedups can be expected when the particle count is very large."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### SPH KERNEL"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The FIRE simulations are meshless finite mass (MFM), not smoothed particle hydrodynamics (SPH) simulations. MARTINI is strictly speaking designed for SPH simulations, but can still provide a good approximation of non-SPH simulations by representing mesh cells as SPH particles assuming their centroids and characteristic smoothing lengths. The MFM scheme in FIRE uses a cubic spline kernel function (although this is not an SPH kernel), so we may as well use MARTINI's `CubicSplineKernel` module for the approximation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sph_kernel = CubicSplineKernel()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### MARTINI"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now set up the configuration:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "M = Martini(\n",
+    "    source=source,\n",
+    "    datacube=datacube,\n",
+    "    beam=beam,\n",
+    "    noise=noise,\n",
+    "    spectral_model=spectral_model,\n",
+    "    sph_kernel=sph_kernel,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Similar to previewing the source, we can make a preview here. Now the extent of the datacube is overlaid with a red box."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "M.preview(\n",
+    "    point_scaling=\"fixed\", vlim=180 * U.km / U.s\n",
+    ")  # can experiment with parameters here"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The main source insertion loop, that is the most computationally demanding step, can be run in parallel if the `multiprocess` package is installed (not to be confused with `multiprocessing`, which is normally included with python!). Edit the cell below to use more than 1 CPU core if this package is installed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbmake": {
+     "post_cell_execute": [
+      "ncpu = 8"
+     ]
+    }
+   },
+   "outputs": [],
+   "source": [
+    "ncpu = 1  # can be >1 if multiprocess package is installed"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we're happy with the preview, we're ready to call the functions to make the actual mock observation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "M.insert_source_in_cube(ncpu=ncpu)\n",
+    "M.add_noise()\n",
+    "M.convolve_beam()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You may notice that the number of pixels in the progress counter differs from the number defined in the DataCube module. This is because convolution with the beam requires some padding, which is ultimately filled with nonsense and discarded.\n",
+    "\n",
+    "To write the results: two output formats are available, depending on preference. Both\n",
+    "formats are self-documenting, via FITS header keywords and HDF5\n",
+    "attributes, respectively. For HDF5 output, the beam image is included\n",
+    "in the same file. (If you do not have h5py installed, comment out the call to `write_hdf5`.)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "M.write_fits(\"fire_martini_demo.fits\")\n",
+    "M.write_beam_fits(\"fire_martini_demo_beam.fits\")\n",
+    "M.write_hdf5(\"fire_martini_demo.hdf5\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Inspect the results (FITS)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's load the FITS file that MARTINI produced and take a quick look."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from astropy.io import fits\n",
+    "from astropy.wcs import WCS\n",
+    "\n",
+    "with fits.open(\"fire_martini_demo.fits\") as f:\n",
+    "    cube_wcs = WCS(f[0].header)\n",
+    "    flux_cube = f[0].data * U.Unit(f[0].header[\"BUNIT\"])\n",
+    "    n_channels = cube_wcs.pixel_shape[cube_wcs.wcs.spec]\n",
+    "    vch = np.array(cube_wcs.sub(axes=[3]).all_pix2world(np.arange(n_channels), 0))[\n",
+    "        0\n",
+    "    ] * U.Unit(cube_wcs.world_axis_units[cube_wcs.wcs.spec])\n",
+    "    vch = vch - source.vsys"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "flux_cube.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's examine one of the velocity channels:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(figsize=(5, 5))\n",
+    "plt.clf()\n",
+    "ax = fig.add_subplot(1, 1, 1, aspect=\"equal\", projection=cube_wcs.celestial)\n",
+    "\n",
+    "# choose units for plotting, not necessarily the units data are stored in:\n",
+    "flux_unit = U.Jy / U.beam\n",
+    "\n",
+    "plt.imshow(flux_cube[27, ...].to_value(flux_unit), cmap=\"Greys\")\n",
+    "plt.xlabel(\"RA\")\n",
+    "plt.ylabel(\"Dec\")\n",
+    "plt.colorbar(label=f\"Flux [{flux_unit}]\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "And do a quick plot of the first three moments:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "# choose plotting units\n",
+    "mom0_unit = U.Jy / U.beam\n",
+    "mom1_unit = U.km / U.s\n",
+    "mom2_unit = U.km / U.s\n",
+    "\n",
+    "rms = np.std(\n",
+    "    flux_cube[:, :20, :20]\n",
+    ")  # noise in a corner patch where there is little signal\n",
+    "clip = np.where(flux_cube > 5 * rms, 1, 0)\n",
+    "np.seterr(all=\"ignore\")\n",
+    "fig = plt.figure(figsize=(16, 5))\n",
+    "plt.clf()\n",
+    "sp1 = fig.add_subplot(1, 3, 1, aspect=\"equal\", projection=cube_wcs.celestial)\n",
+    "sp2 = fig.add_subplot(1, 3, 2, aspect=\"equal\", projection=cube_wcs.celestial)\n",
+    "sp3 = fig.add_subplot(1, 3, 3, aspect=\"equal\", projection=cube_wcs.celestial)\n",
+    "mom0 = np.sum(flux_cube, axis=0)\n",
+    "mask = np.where(mom0 > 0.002 * U.Jy / U.beam, 1, np.nan)\n",
+    "mom1 = np.sum(flux_cube * clip * vch[:, np.newaxis, np.newaxis], axis=0) / mom0\n",
+    "mom2 = np.sqrt(\n",
+    "    np.sum(\n",
+    "        flux_cube\n",
+    "        * clip\n",
+    "        * np.power(vch[:, np.newaxis, np.newaxis] - mom1[np.newaxis], 2),\n",
+    "        axis=0,\n",
+    "    )\n",
+    "    / mom0\n",
+    ")\n",
+    "im1 = sp1.imshow(mom0.to_value(mom0_unit), cmap=\"Greys\")\n",
+    "plt.colorbar(im1, ax=sp1, label=f\"mom0 [{mom0_unit}]\")\n",
+    "im2 = sp2.imshow(\n",
+    "    (mom1 * mask).to_value(mom1_unit),\n",
+    "    cmap=\"RdBu_r\",\n",
+    "    vmin=-np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
+    "    vmax=np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
+    ")\n",
+    "plt.colorbar(im2, ax=sp2, label=f\"mom1 [{mom1_unit}]\")\n",
+    "im3 = sp3.imshow(\n",
+    "    (mom2 * mask).to_value(mom2_unit),\n",
+    "    cmap=\"magma\",\n",
+    ")\n",
+    "plt.colorbar(im3, ax=sp3, label=f\"mom2 [{mom2_unit}]\")\n",
+    "for sp in sp1, sp2, sp3:\n",
+    "    sp.set_xlabel(\"RA\")\n",
+    "    sp.set_ylabel(\"Dec\")\n",
+    "\n",
+    "plt.subplots_adjust(wspace=0.3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This galaxy has a strongly disturbed morphology, but does display clear signs of rotation. The alignment of the disc, such as it is, looks as expected - the inclination is plausibly about 60 degrees, and the position angle of the central part is horizontal in the figure - in this case the automated orientation function has performed well, though it should never be assumed that this will always be the case!\n",
+    "\n",
+    "For complete documentation, more usage examples, and further information, please take a look at the [MARTINI webpage](https://martini.readthedocs.io/en/latest/)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Inspect the results (HDF5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's load the HDF5 that MARTINI produced and take a quick look."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import h5py\n",
+    "\n",
+    "f = h5py.File(\"fire_martini_demo.hdf5\", \"r\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "list(f.keys())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In addition to the mock cube itself (`f[\"FluxCube\"]`), the HDF5 output contains arrays of the same size with the coordinates of the centre of each cell in the cube in RA, Dec and spectral space (`f[\"RA\"]`,`f[\"Dec\"]` and `f[\"channel_mids\"]`). There are also arrays longer by one in each dimension that contain the coordinates of the corners of each cell in the cube (`f[\"RA_vertices\"]`,`f[\"Dec_vertices\"]` and `f[\"channel_vertices\"]`). The latter are convenient for use with the `pcolormesh` plotting function from matplotlib, so we'll read these."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "flux_cube = f[\"FluxCube\"][()] * U.Unit(f[\"FluxCube\"].attrs[\"FluxCubeUnit\"])\n",
+    "ra_vertices = f[\"RA_vertices\"][()] * U.Unit(f[\"RA_vertices\"].attrs[\"Unit\"])\n",
+    "dec_vertices = f[\"Dec_vertices\"][()] * U.Unit(f[\"RA_vertices\"].attrs[\"Unit\"])\n",
+    "spec_vertices = f[\"channel_vertices\"][()] * U.Unit(f[\"channel_vertices\"].attrs[\"Unit\"])\n",
+    "vch = (\n",
+    "    f[\"velocity_channel_mids\"][()] * U.Unit(f[\"velocity_channel_mids\"].attrs[\"Unit\"])\n",
+    "    - source.vsys\n",
+    ")\n",
+    "f.close()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The RA range of our cube straddles the RA=0 boundary, let's shift the `ra_vertices` to a (-180, 180] range for plotting:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ra_vertices = np.where(\n",
+    "    ra_vertices > 180 * U.deg, ra_vertices - 360 * U.deg, ra_vertices\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "flux_cube.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's examine one of the velocity channels:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(figsize=(5, 5))\n",
+    "ax = fig.add_subplot(1, 1, 1, aspect=\"equal\")\n",
+    "\n",
+    "channel_slice = np.s_[:, :, 27]  # a slice of the data cube containing one channel\n",
+    "\n",
+    "# choose units for plotting, not necessarily the units data are stored in:\n",
+    "ra_unit = U.deg\n",
+    "dec_unit = U.deg\n",
+    "flux_unit = U.Jy / U.beam\n",
+    "\n",
+    "plt.pcolormesh(\n",
+    "    ra_vertices[channel_slice].to_value(ra_unit),\n",
+    "    dec_vertices[channel_slice].to_value(dec_unit),\n",
+    "    flux_cube[channel_slice].to_value(flux_unit),\n",
+    "    cmap=\"Greys\",\n",
+    ")\n",
+    "ax.set_xlabel(f\"RA [{ra_unit}]\")\n",
+    "ax.set_ylabel(f\"Dec [{dec_unit}]\")\n",
+    "ax.set_xlim(ax.get_xlim()[::-1])\n",
+    "plt.colorbar(label=f\"Flux [{flux_unit}]\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "And do a quick plot of the first three moments:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "# choose units for plotting, not necessarily the units data are stored in:\n",
+    "ra_unit = U.deg\n",
+    "dec_unit = U.deg\n",
+    "mom0_unit = U.Jy / U.beam\n",
+    "mom1_unit = U.km / U.s\n",
+    "mom2_unit = U.km / U.s\n",
+    "\n",
+    "np.seterr(all=\"ignore\")\n",
+    "fig = plt.figure(figsize=(16, 5))\n",
+    "sp1 = fig.add_subplot(1, 3, 1, aspect=\"equal\")\n",
+    "sp2 = fig.add_subplot(1, 3, 2, aspect=\"equal\")\n",
+    "sp3 = fig.add_subplot(1, 3, 3, aspect=\"equal\")\n",
+    "rms = np.std(\n",
+    "    flux_cube[:16, :16]\n",
+    ")  # noise in a corner patch where there is little signal\n",
+    "clip = np.where(flux_cube > 5 * rms, 1, 0)\n",
+    "mom0 = np.sum(flux_cube, axis=-1)\n",
+    "mask = np.where(mom0 > 0.002 * U.Jy / U.beam, 1, np.nan)\n",
+    "mom1 = np.sum(flux_cube * clip * vch, axis=-1) / mom0\n",
+    "mom2 = np.sqrt(\n",
+    "    np.sum(flux_cube * clip * np.power(vch - mom1[..., np.newaxis], 2), axis=-1) / mom0\n",
+    ")\n",
+    "im1 = sp1.pcolormesh(\n",
+    "    ra_vertices[..., 0].to_value(\n",
+    "        ra_unit\n",
+    "    ),  # pick one channel, coordinates are the same in all of them\n",
+    "    dec_vertices[..., 0].to_value(\n",
+    "        dec_unit\n",
+    "    ),  # pick one channel, coordinates are the same in all of them\n",
+    "    mom0.to_value(mom0_unit),\n",
+    "    cmap=\"Greys\",\n",
+    ")\n",
+    "plt.colorbar(im1, ax=sp1, label=f\"mom0 [{mom0_unit}]\")\n",
+    "im2 = sp2.pcolormesh(\n",
+    "    ra_vertices[..., 0].to_value(\n",
+    "        ra_unit\n",
+    "    ),  # pick one channel, coordinates are the same in all of them\n",
+    "    dec_vertices[..., 0].to_value(\n",
+    "        dec_unit\n",
+    "    ),  # pick one channel, coordinates are the same in all of them\n",
+    "    (mom1 * mask).to_value(mom1_unit),\n",
+    "    cmap=\"RdBu_r\",\n",
+    "    vmin=-np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
+    "    vmax=np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
+    ")\n",
+    "plt.colorbar(im2, ax=sp2, label=f\"mom1 [{mom1_unit}]\")\n",
+    "im3 = sp3.pcolormesh(\n",
+    "    ra_vertices[..., 0].to_value(\n",
+    "        ra_unit\n",
+    "    ),  # pick one channel, coordinates are the same in all of them\n",
+    "    dec_vertices[..., 0].to_value(\n",
+    "        dec_unit\n",
+    "    ),  # pick one channel, coordinates are the same in all of them\n",
+    "    (mom2 * mask).to_value(mom2_unit),\n",
+    "    cmap=\"magma\",\n",
+    "    vmin=0,\n",
+    ")\n",
+    "plt.colorbar(im3, ax=sp3, label=f\"mom2 [{mom2_unit}]\")\n",
+    "for sp in sp1, sp2, sp3:\n",
+    "    sp.set_xlabel(f\"RA [{ra_unit}]\")\n",
+    "    sp.set_ylabel(f\"Dec [{dec_unit}]\")\n",
+    "    sp.set_xlim(sp.get_xlim()[::-1])\n",
+    "plt.subplots_adjust(wspace=0.3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This galaxy has a strongly disturbed morphology, but does display clear signs of rotation. The alignment of the disc, such as it is, looks as expected - the inclination is plausibly about 60 degrees, and the position angle of the central part is horizontal in the figure - in this case the automated orientation function has performed well, though it should never be assumed that this will always be the case!\n",
+    "\n",
+    "For complete documentation, more usage examples, and further information, please take a look at the [MARTINI webpage](https://martini.readthedocs.io/en/latest/)."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/martini_simba.ipynb b/examples/martini_simba.ipynb
index a1b8ff6..de5083e 100644
--- a/examples/martini_simba.ipynb
+++ b/examples/martini_simba.ipynb
@@ -792,7 +792,7 @@
     "plt.colorbar(im1, ax=sp1, label=f\"mom0 [{mom0_unit}]\")\n",
     "im2 = sp2.imshow(\n",
     "    (mom1 * mask).to_value(mom1_unit),\n",
-    "    cmap=\"RdBu\",\n",
+    "    cmap=\"RdBu_r\",\n",
     "    vmin=-np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
     "    vmax=np.nanmax(np.abs(mom1 * mask)).to_value(mom1_unit),\n",
     ")\n",
diff --git a/martini/sources/__init__.py b/martini/sources/__init__.py
index 6d4fd22..79a165f 100644
--- a/martini/sources/__init__.py
+++ b/martini/sources/__init__.py
@@ -11,3 +11,4 @@
 from .magneticum_source import MagneticumSource
 from .simba_source import SimbaSource
 from .tng_source import TNGSource
+from .fire_source import FIRESource
diff --git a/martini/sources/fire_source.py b/martini/sources/fire_source.py
new file mode 100644
index 0000000..f299b94
--- /dev/null
+++ b/martini/sources/fire_source.py
@@ -0,0 +1,193 @@
+import numpy as np
+from ..sph_kernels import _CubicSplineKernel, find_fwhm
+from astropy import units as U, constants as C
+from astropy.coordinates import ICRS
+from .sph_source import SPHSource
+
+
+class FIRESource(SPHSource):
+    """
+    Class abstracting HI sources designed to work with publicly available FIRE
+    snapshot and group data.
+
+    For file access, see https://flathub.flatironinstitute.org/fire
+
+    The authors of the :mod:`gizmo_analysis` package request the folowing: "If you use
+    this package in work that you publish, please cite it, along the lines of: 'This work
+    used GizmoAnalysis (http://ascl.net/2002.015), which first was used in Wetzel et al.
+    2016 (https://ui.adsabs.harvard.edu/abs/2016ApJ...827L..23W).'"
+
+    Parameters
+    ----------
+    simulation_directory : str, optional
+        Base directory containing FIRE simulation output. (Default: ``"."``)
+
+    snapshot_directory : str, optional
+        Directory within ``simulation_directory`` containing snapshots. If ``None``,
+        the default defined by :mod:`gizmo_analysis` is assumed (Default: ``None``)
+
+    snapshot : tuple, optional
+        A 2-tuple specifying the snapshot. The first element is the type of identifier,
+        a string chosen from ``"index"``, ``"redshift"`` or ``"scalefactor"``. The
+        second element is the desired value of the identifier. For example setting
+        ``snapshot=("scalefactor", 0.2)`` will select the closest snapshot to
+        :math:`a=0.2`. (Default: ``("redshift", 0)``)
+
+    host_number : int, optional
+        Galaxy ("host") position in the catalogue, indexed from 1. (Default: ``1``)
+
+    assign_hosts : str, optional
+        Method to compute galaxy centres. Iterative zoom-in based methods are
+        recommended, options include ``"mass"`` (mass-weighted), ``"potential"``
+        (potential-weighted) and ``"massfraction.metals"`` (metallicity-weighted).
+        Other options are ``"track"`` (time-interpolated position) and ``"halo"``
+        (from Rockstar halo catalogue). Setting ``True`` will attempt to find a
+        sensible default by first trying ``"track"`` then ``"mass"``.
+        (Default: ``"mass"``)
+
+    convert_float32 : bool, optional
+        If ``True``, convert floating point values to 32 bit precision to reduce
+        memory usage. (Default: ``False``)
+
+    gizmo_io_verbose : bool, optional
+        If ``True``, allow the gizmo.io module to print progress and diagnostic messages.
+        (Default: ``False``)
+
+    distance : ~astropy.units.Quantity, optional
+        :class:`~astropy.units.Quantity`, with dimensions of length.
+        Source distance, also used to set the velocity offset via Hubble's law.
+        (Default: ``3 * U.Mpc``)
+
+    vpeculiar : ~astropy.units.Quantity, optional
+        :class:`~astropy.units.Quantity`, with dimensions of velocity.
+        Source peculiar velocity, added to the velocity from Hubble's law.
+        (Default: ``0 * U.km * U.s**-1``)
+
+    rotation : dict, optional
+        Must have a single key, which must be one of ``axis_angle``, ``rotmat`` or
+        ``L_coords``. Note that the 'y-z' plane will be the one eventually placed in the
+        plane of the "sky". The corresponding value must be:
+
+        - ``axis_angle`` : 2-tuple, first element one of 'x', 'y', 'z' for the \
+        axis to rotate about, second element a :class:`~astropy.units.Quantity` with \
+        dimensions of angle, indicating the angle to rotate through.
+        - ``rotmat`` : A (3, 3) :class:`~numpy.ndarray` specifying a rotation.
+        - ``L_coords`` : A 2-tuple containing an inclination and an azimuthal \
+        angle (both :class:`~astropy.units.Quantity` instances with dimensions of \
+        angle). The routine will first attempt to identify a preferred plane \
+        based on the angular momenta of the central 1/3 of particles in the \
+        source. This plane will then be rotated to lie in the plane of the \
+        "sky" ('y-z'), rotated by the azimuthal angle about its angular \
+        momentum pole (rotation about 'x'), and inclined (rotation about \
+        'y'). A 3-tuple may be provided instead, in which case the third \
+        value specifies the position angle on the sky (second rotation about 'x'). \
+        The default position angle is 270 degrees.
+
+        (Default: ``np.eye(3)``)
+
+    ra : ~astropy.units.Quantity, optional
+        :class:`~astropy.units.Quantity`, with dimensions of angle.
+        Right ascension for the source centroid. (Default: ``0 * U.deg``)
+
+    dec : ~astropy.units.Quantity, optional
+        :class:`~astropy.units.Quantity`, with dimensions of angle.
+        Declination for the source centroid. (Default: ``0 * U.deg``)
+
+    coordinate_frame : ~astropy.coordinates.builtin_frames.baseradec.BaseRADecFrame, \
+    optional
+        The coordinate frame assumed in converting particle coordinates to RA and Dec, and
+        for transforming coordinates and velocities to the data cube frame. The frame
+        needs to have a well-defined velocity as well as spatial origin. Recommended
+        frames are :class:`~astropy.coordinates.GCRS`, :class:`~astropy.coordinates.ICRS`,
+        :class:`~astropy.coordinates.HCRS`, :class:`~astropy.coordinates.LSRK`,
+        :class:`~astropy.coordinates.LSRD` or :class:`~astropy.coordinates.LSR`. The frame
+        should be passed initialized, e.g. ``ICRS()`` (not just ``ICRS``).
+        (Default: ``astropy.coordinates.ICRS()``)
+    """
+
+    def __init__(
+        self,
+        simulation_directory=".",
+        snapshot_directory=None,
+        snapshot=("redshift", 0),
+        host_number=1,
+        assign_hosts="mass",
+        convert_float32=False,
+        gizmo_io_verbose=False,
+        distance=3.0 * U.Mpc,
+        vpeculiar=0 * U.km / U.s,
+        rotation={"rotmat": np.eye(3)},
+        ra=0.0 * U.deg,
+        dec=0.0 * U.deg,
+        coordinate_frame=ICRS(),
+    ):
+        import gizmo_analysis as gizmo
+
+        gizmo_read_kwargs = dict(
+            species=["gas", "star"],
+            properties=[
+                "position",
+                "velocity",
+                "temperature",
+                "size",
+                "density",
+                "mass",
+                "potential",
+                "massfraction.metals.hydrogen",
+                "hydrogen.neutral.fraction",
+            ],
+            simulation_directory=simulation_directory,
+            snapshot_value_kind=snapshot[0],
+            snapshot_values=snapshot[1],
+            assign_hosts=assign_hosts,
+            host_number=host_number,
+            convert_float32=convert_float32,
+            verbose=gizmo_io_verbose,
+        )
+        if snapshot_directory is not None:
+            gizmo_read_kwargs["snapshot_directory"] = snapshot_directory
+        gizmo_snap = gizmo.io.Read.read_snapshots(
+            **gizmo_read_kwargs,
+        )
+        particles = dict(
+            xyz_g=(
+                gizmo_snap["gas"]["position"]
+                - gizmo_snap.host["position"][host_number - 1]
+            )
+            * gizmo_snap.snapshot["scalefactor"]
+            * U.kpc,
+            vxyz_g=(
+                gizmo_snap["gas"]["velocity"]
+                - gizmo_snap.host["velocity"][host_number - 1]
+            )
+            * U.km
+            * U.s**-1,
+            T_g=gizmo_snap["gas"]["temperature"] * U.K,
+            # see doi:10.1093/mnras/sty1241 Appendix B for molecular partition:
+            mHI_g=np.where(
+                np.logical_and(
+                    gizmo_snap["gas"]["temperature"] * U.K < 300 * U.K,
+                    gizmo_snap["gas"]["density"] * U.Msun * U.kpc**-3
+                    > C.m_p * 10 * U.cm**-3,
+                ),
+                0,
+                gizmo_snap["gas"].prop("mass.hydrogen.neutral"),
+            )
+            * U.Msun,
+            # per A. Wetzel, size / 0.5077 is radius of cpmpact support
+            hsm_g=gizmo_snap["gas"]["size"]
+            / 0.5077
+            * U.kpc
+            * find_fwhm(_CubicSplineKernel().kernel),
+        )
+        super().__init__(
+            **particles,
+            distance=distance,
+            vpeculiar=vpeculiar,
+            rotation=rotation,
+            ra=ra,
+            dec=dec,
+            h=gizmo_snap.info["hubble"],
+            coordinate_frame=coordinate_frame,
+        )
+        return
diff --git a/optional_requirements.txt b/optional_requirements.txt
index 1640edc..a9ae027 100644
--- a/optional_requirements.txt
+++ b/optional_requirements.txt
@@ -12,3 +12,7 @@ g3t
 pytest-cov
 numpydoc
 nbmake
+beautifulsoup4
+gizmo-analysis
+halo-analysis
+utilities-awetzel
diff --git a/pyproject.toml b/pyproject.toml
index 2632f9d..b186e36 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,7 +10,7 @@ authors = [
 ]
 description="Synthetic datacube creation from simulations."
 readme = "README.rst"
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 classifiers = [
     "Programming Language :: Python :: 3",
     "License :: OSI Approved :: GNU General Public License v3 (GPLv3)",
@@ -40,13 +40,18 @@ magneticumsource = ["g3t"]
 simbasource = ["h5py"]
 swiftgalaxysource = [
     "swiftsimio",
-	"swiftgalaxy",
+    "swiftgalaxy",
 ]
 colibresource = [
     "swiftsimio",
     "swiftgalaxy",
     "velociraptor",
 ]
+firesource = [
+    "utilities-awetzel",
+    "gizmo-analysis",
+    "halo-analysis",
+]
 
 [tool.numpydoc_validation]
 checks = [
@@ -80,4 +85,4 @@ exclude = [
        '__init__',
        '__version__',
        '_illustris_tools',
-]
\ No newline at end of file
+]


Test patch:
diff --git a/.github/workflows/lint_and_test.yml b/.github/workflows/lint_and_test.yml
index 3c1bb7f..6bc7a70 100644
--- a/.github/workflows/lint_and_test.yml
+++ b/.github/workflows/lint_and_test.yml
@@ -48,7 +48,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8, 3.9, "3.10", "3.11"]
+        python-version: [3.9, "3.10", "3.11"]
 
     steps:
     - uses: actions/checkout@v4
diff --git a/tests/test_sources.py b/tests/test_sources.py
index b4400fe..f706d90 100644
--- a/tests/test_sources.py
+++ b/tests/test_sources.py
@@ -639,6 +639,51 @@ def test_tng_notebook(self):
                     os.remove(file_to_cleanup)
 
 
+class TestFIRESource:
+    def test_notebook_version(self):
+        """
+        Check that notebook installs current verison of martini.
+        """
+        with open("examples/martini_fire.ipynb") as f:
+            nb_content = f.read()
+        assert (
+            '"!{sys.executable} -m pip install \\"astromartini[firesource]=='
+            + __version__
+            + '\\""'
+            in nb_content
+        )
+
+    @pytest.mark.skipif(
+        not os.path.isfile("examples/m11e_res7100/output/snapshot_600.hdf5"),
+        reason="sample data not locally available",
+    )
+    def test_fire_notebook(self):
+        """
+        Check that the FIRE example notebook runs.
+        """
+        pytest.importorskip(
+            "nbmake", reason="nbmake (optional dependency) not available"
+        )
+        from nbmake.nb_run import NotebookRun
+        from nbmake.pytest_items import NotebookFailedException
+        import pathlib
+
+        nbr = NotebookRun(pathlib.Path("examples/martini_fire.ipynb"), 3600)
+        try:
+            result = nbr.execute()
+            if result.error is not None:
+                raise NotebookFailedException(result)
+        finally:
+            files_to_cleanup = (
+                "examples/fire_martini_demo_beam.fits",
+                "examples/fire_martini_demo.fits",
+                "examples/fire_martini_demo.hdf5",
+            )
+            for file_to_cleanup in files_to_cleanup:
+                if os.path.isfile(file_to_cleanup):
+                    os.remove(file_to_cleanup)
+
+
 class TestMagneticumSource:
     @pytest.mark.xfail
     def test_stuff(self):