{ "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 the input from the [IllustrisTNG](https://ui.adsabs.harvard.edu/abs/2018MNRAS.473.4077P/abstract) simulations. The data are publicly available and hosted at [tng-project.org](https://www.tng-project.org/)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This tutorial can be run either on the [TNG JupyterLab environment](https://www.tng-project.org/data/lab/), or on any system with internet access." ] }, { "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 IllustrisTNG 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.7` 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": {}, "outputs": [], "source": [ "import sys\n", "\n", "!{sys.executable} -m pip install astromartini[tngsource]==2.0.10" ] }, { "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, h5py, six and requests.\n", "\n", "We'll also install matplotlib, used in this notebook:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!{sys.executable} -m pip install matplotlib" ] }, { "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": {}, "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": { "tags": [] }, "source": [ "## TNG Data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `TNGSource` module in MARTINI is designed to run on the IllustrisTNG [JupyterLab](https://www.tng-project.org/data/lab/), or in a standalone mode. If you are running on the TNG JupyterLab then the simulations are stored locally on disk and will be detected and used. Otherwise, the `TNGSource` module will use the TNG [web API](https://www.tng-project.org/data/docs/api/) to download the particles for the galaxy of interest. If running in standalone mode, an API key for the TNG web API is required. You must first [register](https://www.tng-project.org/users/register/) for an account. Once you are registered, your API key can be found [here](https://www.tng-project.org/users/profile/). Enter your API key using the following cell (if you are running this notebook on the TNG JupyterLab, you may leave it blank):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from getpass import getpass\n", "\n", "api_key = getpass(\"TNG web API key:\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can choose a galaxy of interest by selection a [simulation](https://www.tng-project.org/data/docs/background/), a [snapshot](http://www.tng-project.org/data/docs/specifications/#sec1a) and a subhalo ID. One way to search for a subhalo ID is to use [this tool](https://www.tng-project.org/data/search/) - the `ID` column in search results contains subhalo IDs. For example:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "simulation = \"TNG100-1\"\n", "snapshot = 99\n", "subhalo_id = 400547 # a central subhalo with 218 < Vmax < 220, and SFR > 1" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## TNG 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 TNGSource\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, see the [full documentation](https://kyleaoman.github.io/martini/build/html/) for details of all configuration options. A few suggested best-practices specific to TNG are outlined below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SOURCE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that your API key must be provided unless you are using `TNGSource` on the TNG JupyterLab (if you provide it anyways in that case, it will simply be ignored).\n", "\n", "Any downloaded data can be cached using the `cutout_dir` parameter to specify a directory, if a cache for the requested galaxy is found it will be used instead of re-downloading the data. If running on the TNG JupyterLab environment, this parameter is ignored." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "source = TNGSource(\n", " simulation,\n", " snapshot,\n", " subhalo_id,\n", " api_key=api_key,\n", " cutout_dir=\".\",\n", " distance=4 * U.Mpc,\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, so this preview emphasizes the diffuse circumgalactic 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 close to face-on. 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, 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." ] }, { "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": { "tags": [] }, "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=16.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, bmin=30.0 * U.arcsec, bpa=0.0 * U.deg, 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": {}, "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": {}, "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": [ "Since IllustrisTNG uses a moving mesh hydrodynamics solver (Arepo),\n", "there are no formal SPH smoothing lengths and no specified kernel.\n", "However, approximate smoothing lengths can be derived from the cell\n", "volumes and densities, so a reasonable approximation is to use these for imaging. The `TNGSource` module has already computed equivalent SPH smoothing lengths in the correct format for MARTINI, so we just need to choose a smoothing kernel. The `CubicSplineKernel` is perfectly fine here." ] }, { "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()" ] }, { "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()\n", "M.add_noise()\n", "M.convolve_beam()" ] }, { "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!). Simply use `M.insert_source_in_cube(ncpu=N)`, where `N` is the number of processes to run in parallel.\n", "\n", "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(\"tngdemo.fits\", channels=\"velocity\")\n", "M.write_beam_fits(\"tngdemo_beam.fits\", channels=\"velocity\")\n", "M.write_hdf5(\"tngdemo.hdf5\", channels=\"velocity\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inspect the results" ] }, { "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(\"tngdemo.hdf5\", \"r\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "list(f.keys())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FluxCube = f[\"FluxCube\"][()]\n", "vch = f[\"channel_mids\"][()] / 1e3 - source.distance.to(U.Mpc).value * 70 # m/s to km/s\n", "f.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FluxCube.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)\n", "\n", "plt.imshow(FluxCube[:, :, 15].T, cmap=\"Greys\", aspect=1.0, origin=\"lower\")\n", "ax.autoscale(False)\n", "ax.set_xlabel(\"x [px = arcsec/{:.0f}]\".format(datacube.px_size.to(U.arcsec).value))\n", "ax.set_ylabel(\"y [px = arcsec/{:.0f}]\".format(datacube.px_size.to(U.arcsec).value))\n", "plt.colorbar(label=\"Flux [Jy/beam]\");" ] }, { "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", "np.seterr(all=\"ignore\")\n", "fig = plt.figure(figsize=(16, 5))\n", "sp1 = fig.add_subplot(1, 3, 1)\n", "sp2 = fig.add_subplot(1, 3, 2)\n", "sp3 = fig.add_subplot(1, 3, 3)\n", "rms = np.std(FluxCube[:16, :16]) # noise in a corner patch where there is little signal\n", "clip = np.where(FluxCube > 5 * rms, 1, 0)\n", "mom0 = np.sum(FluxCube, axis=-1)\n", "mask = np.where(mom0 > 0.05, 1, np.nan)\n", "mom1 = np.sum(FluxCube * clip * vch, axis=-1) / mom0\n", "mom2 = (\n", " np.sqrt(np.sum(FluxCube * clip * np.power(vch - mom1[..., np.newaxis], 2), axis=-1))\n", " / mom0\n", ")\n", "im1 = sp1.imshow(mom0.T, cmap=\"Greys\", aspect=1.0, origin=\"lower\")\n", "plt.colorbar(im1, ax=sp1, label=\"mom0 [Jy/beam]\")\n", "im2 = sp2.imshow((mom1 * mask).T, cmap=\"RdBu_r\", aspect=1.0, origin=\"lower\")\n", "plt.colorbar(im2, ax=sp2, label=\"mom1 [km/s]\")\n", "im3 = sp3.imshow(\n", " (mom2 * mask).T, cmap=\"magma\", aspect=1.0, origin=\"lower\", vmin=0, vmax=300\n", ")\n", "plt.colorbar(im3, ax=sp3, label=\"mom2 [km/s]\")\n", "for sp in sp1, sp2, sp3:\n", " sp.set_xlabel(\"x [px = arcsec/10]\")\n", " sp.set_ylabel(\"y [px = arcsec/10]\")\n", "plt.subplots_adjust(wspace=0.3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This galaxy clearly has a very nice spiral morphology in HI, has a rotation-dominated velocity field, and is interacting with a companion. The alignment of the disc looks as expected - the inclination looks to be about 60 degrees, and the position angle 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://kyleaoman.github.io/martini)." ] } ], "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.10.0" } }, "nbformat": 4, "nbformat_minor": 4 }