martini module

Create mock observations of HI sources from cosmological hydrodynamical simulations.

Provides Martini, the main class of the package, and a simplified GlobalProfile class for use when only a spectrum (no spatial information) is desired.

class martini.martini.GlobalProfile(*, source: ~martini.sources.sph_source.SPHSource, spectral_model: ~martini.spectral_models._BaseSpectrum, n_channels: int, channel_width: ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("km / s")] | ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("Hz")], spectral_centre: ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("km / s")] | ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("Hz")] = <Quantity 0. km / s>, quiet: bool = False, channels: None = None)[source]

Bases: _BaseMartini

Just produce a global spectrum of an HI source.

A simplified version of the main Martini class for those times when you just want a spectrum.

Sometimes only the spatially integrated spectrum of a source is wanted, which makes many parts of the standard MARTINI process unnecessary. This class streamlines the process of producing a spatially-integrated spectrum, or “global profile”, with some specific assumptions:

The GlobalProfile class does not assume any spatial aperture, instead every particle in the source contributes to the spectrum (unless it falls outside of the spectral bandwidth).

The positions of particles are still used to calculate the line-of-sight vector and the velocity along this direction.

There is therefore no need or way to specify a beam or SPH kernel as with the main Martini class. It is also not possible to use MARTINI’s noise modules with this class. If these restrictions are found to be too limiting (for example, if the spectrum within a spatial mask defined by a signal-to-noise or other cut is desired), the best course of action is to produce a spatially-resolved mock observation and derive the spectrum from those data as would be done with “real” observations. This class is mainly intended to efficiently provide a “quick look” at the spectrum, or a reference “ideal” spectrum.

Parameters:

source (SPHSource) – An instance of a class derived from SPHSource. A description of the HI emitting object, including position, geometry and an interface to the simulation data (SPH particle masses, positions, etc.). See sub-module documentation.
spectral_model (_BaseSpectrum) – An instance of a class derived from _BaseSpectrum. A description of the HI line produced by a particle of given properties. A Dirac-delta spectrum, and both fixed-width and temperature-dependent Gaussian line models are provided; implementing other models is straightforward. See sub-module documentation.
n_channels (int) – Number of channels along the spectral axis.
channel_width (Quantity) – Quantity, with dimensions of velocity or frequency. Step size along the spectral axis. Can be provided as a velocity or a frequency.
spectral_centre (Quantity, optional) – Quantity with dimensions of velocity or frequency. Velocity (or frequency) of the centre along the spectral axis. Defaults to 0 km/s, but should normally be set to the systemic velocity of the source.
quiet (bool, optional) – If True, suppress output to stdout.
channels (str, deprecated) – Deprecated, channels and their units now fixed at DataCube initialization.

Examples

# ------make a toy galaxy----------
N = 500
phi = np.random.rand(N) * 2 * np.pi
r = []
for L in np.random.rand(N):

    def f(r):
        return L - 0.5 * (2 - np.exp(-r) * (np.power(r, 2) + 2 * r + 2))

    r.append(fsolve(f, 1.0)[0])
r = np.array(r)
# exponential disk
r *= 3 / np.sort(r)[N // 2]
z = -np.log(np.random.rand(N))
# exponential scale height
z *= 0.5 / np.sort(z)[N // 2] * np.sign(np.random.rand(N) - 0.5)
x = r * np.cos(phi)
y = r * np.sin(phi)
xyz_g = np.vstack((x, y, z)) * U.kpc
# linear rotation curve
vphi = 100 * r / 6.0
vx = -vphi * np.sin(phi)
vy = vphi * np.cos(phi)
# small pure random z velocities
vz = (np.random.rand(N) * 2.0 - 1.0) * 5
vxyz_g = np.vstack((vx, vy, vz)) * U.km * U.s**-1
T_g = np.ones(N) * 8e3 * U.K
mHI_g = np.ones(N) / N * 5.0e9 * U.Msun
# ~mean interparticle spacing smoothing
hsm_g = np.ones(N) * 4 / np.sqrt(N) * U.kpc
# ---------------------------------

source = SPHSource(
    distance=3.0 * U.Mpc,
    rotation={"L_coords": (60.0 * U.deg, 0.0 * U.deg)},
    ra=0.0 * U.deg,
    dec=0.0 * U.deg,
    h=0.7,
    T_g=T_g,
    mHI_g=mHI_g,
    xyz_g=xyz_g,
    vxyz_g=vxyz_g,
    hsm_g=hsm_g,
)

spectral_model = GaussianSpectrum(sigma=7 * U.km * U.s**-1)

M = GlobalProfile(
    source=source,
    spectral_model=spectral_model,
    n_channels=32,
    channel_width=10 * U.km * U.s**-1,
    spectral_centre=source.vsys,
)

# spectrum and channel edges and centres can be accessed as:
M.spectrum  # spectrum is evaluated at first access if not already done
M.channel_edges
M.channel_mids

# the spectrum can be explicitly evaluated with:
M.insert_source_in_spectrum()

property channel_edges: Annotated[Quantity, Unit('Hz')] | Annotated[Quantity, Unit('km / s')]

The edges of the channels for the spectrum.

Returns:: Quantity with dimensions of frequency or velocity. Edges of the channels with units depending on GlobalProfile’s native channel spacing.
Return type:: Quantity

property channel_mids: Annotated[Quantity, Unit('Hz')] | Annotated[Quantity, Unit('km / s')]

The centres of the channels for the spectrum.

Returns:: Quantity with dimensions of frequency or velocity. Edges of the channels with units depending on GlobalProfile’s native channel spacing.
Return type:: Quantity

property channel_width: Annotated[Quantity, Unit('Hz')] | Annotated[Quantity, Unit('km / s')]

The width of the channels for the spectrum.

Returns:: Quantity with dimentions of frequency or velocity. Width of the channels with units depending on GlobalProfile’s channels argument.
Return type:: Quantity

property frequency_channel_edges: Annotated[Quantity, Unit('Hz')]

The edges of the frequency channels for the spectrum.

Returns:: Quantity with dimensions of frequency. Edges of the channels with frequency units.
Return type:: Quantity

property frequency_channel_mids: Annotated[Quantity, Unit('Hz')]

The centres of the frequency channels for the spectrum.

Returns:: Quantity with dimensions of frequency. Edges of the channels with frequency units.
Return type:: Quantity

init_spectra() → None

Explicitly trigger evaluation of particle spectra.

Triggered when insert_source_in_cube() is called if not called explicitly before this.

insert_source_in_spectrum() → None[source]

Populate the DataCube with flux from source particles.

For the GlobalProfile class we assume that every particle in the source contributes (if it falls in the spectral bandwidth) regardless of position on the sky. The line-of-sight vector still depends on the particle positions, so the direction to the individual particles is still taken into account.

plot_spectrum(fig: int = 1, title: str = '', channels: str = 'velocity', show_vsys: bool = True, save: str | None = None) → Figure[source]

Produce a figure showing the spectrum.

Parameters:

fig (int, optional) – Number of the figure in matplotlib, it will be created as plt.figure(fig).
title (str, optional) – A title for the figure can be provided.
channels (str, optional) – The type of spectral axis for the plot, either "velocity" or "frequency".
show_vsys (bool, optional) – If True, draw a vertical line at the source systemic velocity.
save (str, optional) – If provided, the figure is saved using plt.savefig(save). A .png or .pdf suffix is recommended.

Returns:

The spectrum Figure.

Return type:

matplotlib.figure.Figure

preview(max_points: int = 5000, fig: int = 1, lim: str | Annotated[Quantity, Unit('deg')] | None = None, vlim: str | Annotated[Quantity, Unit('km / s')] | None = None, point_scaling: str = 'auto', title: str = '', save: str | None = None) → Figure

Produce a figure showing the source particle coordinates and velocities.

The data cube region is also outlined.

Makes a 3-panel figure showing the projection of the source as it will appear in the mock observation. The first panel shows the particles in the y-z plane, coloured by the x-component of velocity (MARTINI projects the source along the x-axis). The second and third panels are position-velocity diagrams showing the x-component of velocity against the y and z coordinates, respectively. The red boxes drawn on the panels show the data cube extent in position and velocity at the distance and systemic velocity of the source.

Parameters:

max_points (int, optional) – Maximum number of points to draw per panel, the particles will be randomly subsampled if the source has more.
fig (int, optional) – Number of the figure in matplotlib, it will be created as plt.figure(fig).
lim (Quantity, optional) – Quantity with dimensions of angle. The coordinate axes extend from -lim to lim. If unspecified, the maximum angular separation between the source centre defined by its ra and dec and the particle sky positions is used. Whether specified or not, the axes are expanded if necessary to contain the extent of the data cube. Alternatively, the string "datacube" can be passed. In this case the axis limits are set by the extent of the data cube.
vlim (Quantity, optional) – Quantity with dimensions of speed. The velocity axes and colour bar extend from -vlim to vlim. If unspecified, the maximum absolute velocity of particles in the source is used. Whether specified or not, the axes are expanded if necessary to contain the extent of the data cube. Alternatively, "datacube" can be passed. In this case the axis limits are set by the extent of the data cube.
point_scaling (str, optional) – By default points are scaled in size and transparency according to their HI mass and the smoothing length (loosely proportional to their surface densities, but with different scaling to achieve a visually useful plot). For some sources the automatic scaling may not give useful results, using point_scaling="fixed" will plot points of constant size without opacity.
title (str, optional) – A title for the figure can be provided.
save (str, optional) – If provided, the figure is saved using plt.savefig(save). A .png or .pdf suffix is recommended.

Returns:

The preview Figure.

Return type:

matplotlib.figure.Figure

reset() → None[source]: Re-initialize the spectrum with zero-values.

property spectrum: Annotated[Quantity, Unit('Jy')]

The spectrum of the source with spatial information integrated out.

The spectrum is calculated when this property is accessed if this hasn’t already been done.

Returns:: Quantity with dimensions of flux density. Spatially-integrated spectrum of the source.
Return type:: Quantity

property velocity_channel_edges: Annotated[Quantity, Unit('km / s')]

The edges of the channels for the spectrum in velocity units.

Returns:: Quantity with dimensions of velocity. Edges of the channels with velocity units.
Return type:: Quantity

property velocity_channel_mids: Annotated[Quantity, Unit('km / s')]

The centres of the channels for the spectrum.

Returns:: Quantity with dimensions of velocity. Edges of the channels with velocity units.
Return type:: Quantity

class martini.martini.Martini(*, source: SPHSource, datacube: DataCube, beam: _BaseBeam | None = None, noise: _BaseNoise | None = None, sph_kernel: _BaseSPHKernel, spectral_model: _BaseSpectrum, quiet: bool = False)[source]

Bases: _BaseMartini

Creates synthetic HI data cubes from simulation data.

Usual use of martini involves first creating instances of classes from each of the required and optional sub-modules, then creating a Martini with these instances as arguments. The object can then be used to create synthetic observations, usually by calling insert_source_in_cube(), (optionally) add_noise(), (optionally) convolve_beam() and write_fits() in order.

Parameters:

source (SPHSource) – An instance of a class derived from SPHSource. A description of the HI emitting object, including position, geometry and an interface to the simulation data (SPH particle masses, positions, etc.). See sub-module documentation.
datacube (DataCube) – A DataCube instance. A description of the datacube to create, including pixels, channels, sky position. See sub-module documentation.
beam (_BaseBeam, optional) – An instance of a class derived from ~martini.beams._BaseBeam. A description of the beam for the simulated telescope. Given a description, either mathematical or as an image, the creation of a custom beam is straightforward. See sub-module documentation.
noise (_BaseNoise, optional) – An instance of a class derived from _BaseNoise. A description of the simulated noise. A simple Gaussian noise model is provided; implementation of other noise models is straightforward. See sub-module documentation.
sph_kernel (_BaseSPHKernel) – An instance of a class derived from _BaseSPHKernel. A description of the SPH smoothing kernel. Check simulation documentation for the kernel used in a particular simulation, and SPH kernel sub-module documentation for guidance.
spectral_model (_BaseSpectrum) – An instance of a class derived from _BaseSpectrum. A description of the HI line produced by a particle of given properties. A Dirac-delta spectrum, and both fixed-width and temperature-dependent Gaussian line models are provided; implementing other models is straightforward. See sub-module documentation.
quiet (bool, optional) – If True, suppress output to stdout.

Examples

More detailed examples can be found in the examples directory in the github distribution of the package.

The following example illustrates basic use of MARTINI, using a (very!) crude model of a gas disk. This example can be run by doing ‘from martini import demo; demo()’:

# ------make a toy galaxy----------
N = 500
phi = np.random.rand(N) * 2 * np.pi
r = []
for L in np.random.rand(N):

    def f(r):
        return L - 0.5 * (2 - np.exp(-r) * (np.power(r, 2) + 2 * r + 2))

r.append(fsolve(f, 1.0)[0])
r = np.array(r)
# exponential disk
r *= 3 / np.sort(r)[N // 2]
z = -np.log(np.random.rand(N))
# exponential scale height
z *= 0.5 / np.sort(z)[N // 2] * np.sign(np.random.rand(N) - 0.5)
x = r * np.cos(phi)
y = r * np.sin(phi)
xyz_g = np.vstack((x, y, z)) * U.kpc
# linear rotation curve
vphi = 100 * r / 6.0
vx = -vphi * np.sin(phi)
vy = vphi * np.cos(phi)
# small pure random z velocities
vz = (np.random.rand(N) * 2.0 - 1.0) * 5
vxyz_g = np.vstack((vx, vy, vz)) * U.km * U.s**-1
T_g = np.ones(N) * 8e3 * U.K
mHI_g = np.ones(N) / N * 5.0e9 * U.Msun
# ~mean interparticle spacing smoothing
hsm_g = np.ones(N) * 4 / np.sqrt(N) * U.kpc
# ---------------------------------

source = SPHSource(
    distance=3.0 * U.Mpc,
    rotation={"L_coords": (60.0 * U.deg, 0.0 * U.deg)},
    ra=0.0 * U.deg,
    dec=0.0 * U.deg,
    h=0.7,
    T_g=T_g,
    mHI_g=mHI_g,
    xyz_g=xyz_g,
    vxyz_g=vxyz_g,
    hsm_g=hsm_g,
)

datacube = DataCube(
    n_px_x=128,
    n_px_y=128,
    n_channels=32,
    px_size=10.0 * U.arcsec,
    channel_width=10.0 * U.km * U.s**-1,
    spectral_centre=source.vsys,
)

beam = GaussianBeam(
    bmaj=30.0 * U.arcsec, bmin=30.0 * U.arcsec, bpa=0.0 * U.deg, truncate=4.0
)

noise = GaussianNoise(rms=3.0e-5 * U.Jy * U.beam**-1)

spectral_model = GaussianSpectrum(sigma=7 * U.km * U.s**-1)

sph_kernel = CubicSplineKernel()

M = Martini(
    source=source,
    datacube=datacube,
    beam=beam,
    noise=noise,
    spectral_model=spectral_model,
    sph_kernel=sph_kernel,
)

M.insert_source_in_cube()
M.add_noise()
M.convolve_beam()
M.write_beam_fits(beamfile)
M.write_fits(cubefile)
print(f"Wrote demo fits output to {cubefile}, and beam image to {beamfile}.")
try:
    M.write_hdf5(hdf5file)
except ModuleNotFoundError:
    print("h5py package not present, skipping hdf5 output demo.")
else:
    print(f"Wrote demo hdf5 output to {hdf5file}.")

add_noise() → None[source]: Insert noise into the data cube.

convolve_beam() → None[source]: Convolve the beam and data cube.

property datacube: DataCube

The DataCube object for this mock observation.

Returns:: The DataCube contained by this Martini instance.
Return type:: DataCube

init_spectra() → None

Explicitly trigger evaluation of particle spectra.

Triggered when insert_source_in_cube() is called if not called explicitly before this.

insert_source_in_cube(skip_validation: bool = False, progressbar: bool | None = None, ncpu: int = 1) → None[source]

Populate the DataCube with flux from the particles in the source.

Parameters:

skip_validation (bool, optional) – SPH kernel interpolation onto the DataCube is approximated for increased speed. For some combinations of pixel size, distance and SPH smoothing length, the approximation may break down. The kernel class will check whether this will occur and raise a RuntimeError if so. This validation can be skipped (at the cost of accuracy!) by setting this parameter True.
progressbar (bool, optional) – A progress bar is shown by default. Progress bars work, with perhaps some visual glitches, in parallel. If Martini was initialised with quiet set to True, progress bars are switched off unless explicitly turned on.
ncpu (int) – Number of processes to use in main source insertion loop. Using more than one cpu requires the multiprocess module (n.b. not the same as multiprocessing).

preview(max_points: int = 5000, fig: int = 1, lim: str | Annotated[Quantity, Unit('deg')] | None = None, vlim: str | Annotated[Quantity, Unit('km / s')] | None = None, point_scaling: str = 'auto', title: str = '', save: str | None = None) → Figure

Produce a figure showing the source particle coordinates and velocities.

The data cube region is also outlined.

Makes a 3-panel figure showing the projection of the source as it will appear in the mock observation. The first panel shows the particles in the y-z plane, coloured by the x-component of velocity (MARTINI projects the source along the x-axis). The second and third panels are position-velocity diagrams showing the x-component of velocity against the y and z coordinates, respectively. The red boxes drawn on the panels show the data cube extent in position and velocity at the distance and systemic velocity of the source.

Parameters:

max_points (int, optional) – Maximum number of points to draw per panel, the particles will be randomly subsampled if the source has more.
fig (int, optional) – Number of the figure in matplotlib, it will be created as plt.figure(fig).
lim (Quantity, optional) – Quantity with dimensions of angle. The coordinate axes extend from -lim to lim. If unspecified, the maximum angular separation between the source centre defined by its ra and dec and the particle sky positions is used. Whether specified or not, the axes are expanded if necessary to contain the extent of the data cube. Alternatively, the string "datacube" can be passed. In this case the axis limits are set by the extent of the data cube.
vlim (Quantity, optional) – Quantity with dimensions of speed. The velocity axes and colour bar extend from -vlim to vlim. If unspecified, the maximum absolute velocity of particles in the source is used. Whether specified or not, the axes are expanded if necessary to contain the extent of the data cube. Alternatively, "datacube" can be passed. In this case the axis limits are set by the extent of the data cube.
point_scaling (str, optional) – By default points are scaled in size and transparency according to their HI mass and the smoothing length (loosely proportional to their surface densities, but with different scaling to achieve a visually useful plot). For some sources the automatic scaling may not give useful results, using point_scaling="fixed" will plot points of constant size without opacity.
title (str, optional) – A title for the figure can be provided.
save (str, optional) – If provided, the figure is saved using plt.savefig(save). A .png or .pdf suffix is recommended.

Returns:

The preview Figure.

Return type:

matplotlib.figure.Figure

reset() → None: Re-initialize the DataCube with zero-values.

write_beam_fits(filename: str, overwrite: bool = True, channels: None = None) → None[source]

Output the beam to a FITS-format file.

The beam is written to file, with pixel sizes, coordinate system, etc. similar to those used for the DataCube.

Parameters:

filename (str) – Name of the file to write. ".fits" will be appended if not already present.
overwrite (bool, optional) – Whether to allow overwriting existing files.
channels (str, deprecated) – Deprecated, channels and their units now fixed at DataCube initialization.

Raises:

ValueError – If Martini was initialized without a beam.

write_fits(filename: str, overwrite: bool = True, obj_name: str = 'MOCK', channels: None = None, dtype: str | dtype | None = None) → None[source]

Output the data cube to a FITS-format file.

Parameters:

filename (str) – Name of the file to write. '.fits' will be appended if not already present in the filename. '.fits.gz' etc. will be recognized.
overwrite (bool, optional) – Whether to allow overwriting existing files.
obj_name (str) – Name to write in the OBJECT FITS header field (max 16 characters).
channels (str, deprecated) – Deprecated, channels and their units now fixed at DataCube initialization.
dtype (str or dtype) – Typecode or data-type to which the array is cast. Should be supported by fits. Default to not do data type conversion.

write_hdf5(filename: str, overwrite: bool = True, memmap: bool = False, compact: bool = False, channels: None = None) → h5py.File | None[source]

Output the data cube and beam to a HDF5-format file.

Requires the h5py package.

Parameters:

filename (str) – Name of the file to write. '.hdf5' will be appended if not already present.
overwrite (bool, optional) – Whether to allow overwriting existing files.
memmap (bool, optional) – If True, create a file-like object in memory and return it instead of writing file to disk.
compact (bool, optional) – If True, omit pixel coordinate arrays to save disk space. In this case pixel coordinates can still be reconstructed from FITS-style keywords stored in the FluxCube attributes.
channels (str, deprecated) – Deprecated, channels and their units now fixed at DataCube initialization.