martini package

Mock APERTIF-like Radio Telescope Interferometry of the Neutral ISM.

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 then 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.

class martini.DataCube(*, n_px_x: int, n_px_y: int, n_channels: int, px_size: ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("arcsec")], channel_width: ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("km / s")], spectral_centre: ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("km / s")] = <Quantity 0. km / s>, ra: ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("deg")] = <Quantity 0. deg>, dec: ~typing.Annotated[~astropy.units.quantity.Quantity, Unit("deg")] = <Quantity 0. deg>, stokes_axis: bool = False, coordinate_frame: BaseRADecFrame = <ICRS Frame>, specsys: str = 'icrs', velocity_centre: None = None)[source]

Bases: object

Handles creation and management of the data cube itself.

Basic usage simply involves initializing with the parameters listed below. More advanced usage might arise if designing custom classes for other sub- modules, especially beams. To initialize a DataCube from a saved state, see load_state().

Parameters:

n_px_x (int) – Pixel count along the x (RA) axis. Even integers strongly preferred.
n_px_y (int) – Pixel count along the y (Dec) axis. Even integers strongly preferred.
n_channels (int) – Number of channels along the spectral axis.
px_size (Quantity) – Quantity, with dimensions of angle. Angular scale of one pixel.
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, normally set this to the systemic velocity of the source.
ra (Quantity, optional) – Quantity, with dimensions of angle. Right ascension of the cube centroid. Defaults to 0 degrees.
dec (Quantity, optional) – Quantity with dimensions of angle. Declination of the cube centroid. Defaults to 0 degrees.
stokes_axis (bool, optional) – Whether the datacube should be initialized with a Stokes’ axis.
coordinate_frame (BaseRADecFrame, optional) – The coordinate frame of the World Coordinate System (WCS) associated with the data cube. Recommended frames are GCRS, ICRS, HCRS, LSRK, LSRD or LSR. The frame should be passed initialized, e.g. ICRS() (not just ICRS).
specsys (str, optional) – The spectral reference frame (standard of rest) of the World Coordinate System (WCS) associated with the data cube. Some common options include "gcrs", "icrs", "hcrs", "lsrk", "lsrd", "lsr". For a complete list, use astropy.coordinates.frame_transform_graph.get_names().
velocity_centre (Quantity, deprecated) – Deprecated, use spectral centre instead.

Notes

The channel_width argument defines the channel spacing in either frequency or velocity units. If provided with units of frequency, the data cube channels will be evenly spaced in frequency. Conversely, if provided with units of velocity, the channels will be evenly spaced in velocity. The spectral_centre is related and can also have units of frequency or velocity, but will be converted to have the same units as the channel_width in order to define the channels.

add_pad(pad: tuple[int, int]) → None[source]

Resize the cube to add a padding region in the spatial direction.

Accurate convolution with a beam requires a cube padded according to the size of the beam kernel (its representation sampled on a grid with the same spacing). The beam class is required to handle defining the size of pad required.

Parameters:: pad (tuple) – 2-tuple (or other sequence) containing the number of pixels to add in the x (RA) and y (Dec) directions.

See also

martini.datacube.DataCube.drop_pad

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

The edges of the channels from the coordinate system.

Returns:: Quantity with dimensions of frequency or velocity containing the channel edges.
Return type:: Quantity

property channel_maps: Iterator[Quantity]

An iterator over the channel maps.

An alias for spatial_slices().

Returns:: The iterator over the spatial ‘slices’ of the cube.
Return type:: iter

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

The centres of the channels from the coordinate system.

Returns:: Quantity with dimensions of frequency or velocity containing the channel centres.
Return type:: Quantity

copy() → Self[source]

Produce a copy of the DataCube.

May be especially useful to create multiple datacubes with differing intermediate steps.

Returns:: Copy of the DataCube object.
Return type:: DataCube

drop_pad() → None[source]

Remove the padding added using add_pad().

After convolution, the pad region contains meaningless information and can be discarded.

See also

martini.datacube.DataCube.add_pad

freq_channels() → None[source]: Issue a warning then do nothing (deprecated).

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

The edges of the channels from the coordinate system in frequency units.

Returns:: Quantity with dimensions of frequency containing the channel edges.
Return type:: Quantity

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

The centres of the channels from the coordinate system in frequency units.

Returns:: Quantity with dimensions of frequency containing the channel centres.
Return type:: Quantity

classmethod from_wcs(input_wcs: WCS, specsys: str | None = None) → Self[source]

Create a DataCube from a World Coordinate System (WCS).

The WCS can be obtained, for instance, from a FITS header.

To create a MARTINI data cube with pixels and channels exactly matching an observed data cube (stored in a FITS file) would be a bit tedious using the usual __init__() since the number of pixels, channel spacing and so on must be specified by hand. This function instead constructs a DataCube from a WCS instance, that can be easily obtained from a FITS file or header (see example below). The resulting data cube has exactly the same dimensions (pixels and channels) and World Coordinate System (WCS) as the input WCS.

Parameters:

input_wcs (WCS) – The WCS instance to use as the basis for the DataCube.
specsys (str, optional) – The spectral reference frame (standard of rest) of the World Coordinate System (WCS) associated with the data cube, selected from the list "gcrs", "icrs", "hcrs", "lsrk", "lsrd", "lsr".

See also

martini.datacube.DataCube

Notes

MARTINI’s data cubes have a fixed axis ordering: first the RA axis, then the Dec axis, then the spectral axis, and finally the Stokes’ axis (if present). A DataCube created with from_wcs() may therefore have its axes re-ordered (tranposed) relative to the input_wcs. The FITS files output by MARTINI have the same axis ordering, and may therefore also be transposed relative to a data cube used to construct the WCS for the input_wcs argument.

Examples

It is easy to initialize a WCS from a FITS-format header. For example, given a FITS file my_cube.fits, setting up a DataCube with matching World Coordinate System (WCS) looks like:

from astropy import wcs
from astropy.io import fits
from martini.datacube import DataCube

with fits.open("my_cube.fits") as fitsfile:
    fits_hdr = fitsfile[0].header  # header of the main HDU
fits_wcs = wcs.WCS(fits_hdr)
datacube = DataCube.from_wcs(fits_wcs)

classmethod load_state(filename: str) → Self[source]

Initialize a DataCube from a state saved.

Note that h5py must be installed for use. Note that ONLY the DataCube state is restored, other modules and their configurations are not affected.

Parameters:: filename (str) – File to open.
Returns:: A suitably initialized DataCube object.
Return type:: DataCube

See also

martini.datacube.DataCube.save_state

save_state(filename: str, overwrite: bool = False) → None[source]

Write the DataCube state to file for re-use.

The state can be re-initialized (see load_state()). Note that h5py must be installed for use. NOT for outputting mock observations, for this see write_fits() and write_hdf5().

Parameters:

filename (str) – File to write.
overwrite (bool) – Whether to allow overwriting existing files. (default: False).

See also

martini.datacube.DataCube.load_state

property spatial_slices: Iterator[Quantity]

An iterator over the spatial ‘slices’ of the cube.

Returns:: The iterator over the spatial ‘slices’ of the cube.
Return type:: iter

property spectra: Iterator[Quantity]

An iterator over the spectra (one in each spatial pixel).

Returns:: The iterator over the spectra making up the cube.
Return type:: iter

property units: tuple[Annotated[Quantity, Unit('deg')], Annotated[Quantity, Unit('deg')], Annotated[Quantity, Unit('Hz')] | Annotated[Quantity, Unit('m / s')]] | tuple[Annotated[Quantity, Unit('deg')], Annotated[Quantity, Unit('deg')], Annotated[Quantity, Unit('Hz')] | Annotated[Quantity, Unit('m / s')], Annotated[Quantity, Unit(dimensionless)]]

The units of the DataCube’s World Coordinate System (WCS).

Returns:: A tuple containing Unit instances describing the WCS units.
Return type:: tuple

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

The edges of the channels from the coordinate system in velocity units.

Returns:: Quantity with dimensions of velocity containing the channel edges.
Return type:: Quantity

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

The centres of the channels from the coordinate system in velocity units.

Returns:: Quantity with dimensions of velocity containing the channel centres.
Return type:: Quantity

velocity_channels() → None[source]: Issue a warning then do nothing (deprecated).

property wcs: WCS

The DataCube’s World Coordinate System (WCS).

Returns:: The WCS instance that describes the DataCube’s World Coordinate System (WCS).
Return type:: WCS

class 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.L_coords(incl: Quantity, Unit("deg")]=<Quantity 0. deg>, az_rot: Quantity, Unit("deg")]=<Quantity 0. deg>, pa: Quantity, Unit("deg")]=<Quantity 270. deg>)[source]

Bases: NamedTuple

Provide an unambiguous way to specify a source orientation based on angular momentum.

The orientation is defined as follows. First the angular momentum vector of the central 1/3 of particles weighted by HI mass is calculated. The angular momentum vector is then rotated to point along the x axis (the line of sight). The disc is now face-on. The L_coords arguments are then applied. First the source can be rotated around its angular momentum vector by an angle az_rot. Then it can be rotated to incline it to the line of sight by an angle incl. Finally it can be rotated in the plane of the sky by an angle pa. All rotations are right-handed.

The initialisation arguments are:

incl (Quantity, optional): The inclination with units of angle, defaults to 0 degrees (face-on).

az_rot (Quantity, optional): The rotation about the angular momentum axis with units of angle, defaults to 0 degrees.

pa (Quantity, optional): The position angle on the sky with units of angle, defaults to 270 degrees.

az_rot: Annotated[Quantity, Unit('deg')]: Alias for field number 1

count(value, /): Return number of occurrences of value.

incl: Annotated[Quantity, Unit('deg')]: Alias for field number 0

index(value, start=0, stop=9223372036854775807, /)

Return first index of value.

Raises ValueError if the value is not present.

pa: Annotated[Quantity, Unit('deg')]: Alias for field number 2

class 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.

martini.demo(cubefile: str = 'testcube.fits', beamfile: str = 'testbeam.fits', hdf5file: str = 'testcube.hdf5') → None[source]

Demonstrate basic usage of MARTINI.

Creates a (very!) crude toy model of a galaxy with a linearly rising rotation curve, exponential disk profile, exponential vertical structure. A basic configuration of MARTINI is initialized and used to create and output a datacube and an image of the beam.

Parameters:

cubefile (str) – File name to write demonstration datacube.
beamfile (str) – File name to write demonstration beam.
hdf5file (str) – File name to write demonstration datacube in hdf5 format (if h5py is available).

martini.demo_source(N: int = 500) → SPHSource[source]

Create a simple toy model of a galaxy.

Parameters:: N (int) – Number of particles to generate in source.
Returns:: An initialized MARTINI source module containing a toy model of a galaxy.
Return type:: SPHSource

martini package

Submodules