martini module
- class martini.martini.GlobalProfile(source=None, spectral_model=None, n_channels=64, channel_width=<Quantity 4. km / s>, velocity_centre=<Quantity 0. km / s>, channels='frequency', quiet=False)[source]
Bases:
_BaseMartini
A simplified version of the main
Martini
class to just produce 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, optional) – Number of channels along the spectral axis. (Default:
64
)channel_width (Quantity, optional) –
Quantity
, with dimensions of velocity or frequency. Step size along the spectral axis. Can be provided as a velocity or a frequency. (Default:4 * U.km / U.s
)velocity_centre (Quantity, optional) –
Quantity
with dimensions of velocity or frequency. Velocity (or frequency) of the centre along the spectral axis. (Default:0 * U.km / U.s
)channels (str, optional) – Type of units (
"velocity"
or"frequency"
) used along the spectral axis. (Default:"frequency"
)quiet (bool) – If
True
, suppress output to stdout. (Default:False
)
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, velocity_centre=source.vsys, channels="frequency", ) # 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()
See also
- property channel_edges
The edges of the channels for the spectrum.
- Returns:
out –
Quantity
with dimensions of frequency or velocity. Edges of the channels with units depending onGlobalProfile
’schannels
argument.- Return type:
See also
- property channel_mids
The centres of the channels for the spectrum.
- Returns:
out –
Quantity
with dimensions of frequency or velocity. Edges of the channels with units depending onGlobalProfile
’schannels
argument.- Return type:
See also
- property channel_width
The width of the channels for the spectrum.
- Returns:
out –
Quantity
with dimentions of frequency or velocity. Width of the channels with units depending onGlobalProfile
’schannels
argument.- Return type:
- insert_source_in_spectrum()[source]
Populates the
DataCube
with flux from the particles in the source.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=1, title='', show_vsys=True, save=None)[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)
. (Default:1
)title (str, optional) – A title for the figure can be provided. (Default:
""
)show_vsys (bool, optional) – If
True
, draw a vertical line at the source systemic velocity. (Default:True
)save (str, optional) – If provided, the figure is saved using
plt.savefig(save)
. A .png or .pdf suffix is recommended. (Default:None
)
- Returns:
out – The spectrum
Figure
.- Return type:
Figure
- preview(max_points=5000, fig=1, lim=None, vlim=None, point_scaling='auto', title='', save=None)
Produce a figure showing the source particle coordinates and velocities, and the data cube region.
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. (Default:
1000
)fig (int, optional) – Number of the figure in matplotlib, it will be created as
plt.figure(fig)
. (Default:1
)lim (Quantity, optional) –
Quantity
with dimensions of length. The coordinate axes extend from-lim
tolim
. If unspecified, the maximum absolute coordinate 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, the string"datacube"
can be passed. In this case the axis limits are set by the extent of the data cube. (Default:None
)vlim (Quantity, optional) –
Quantity
with dimensions of speed. The velocity axes and colour bar extend from-vlim
tovlim
. 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. (Default:None
)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. (Default:"auto"
)title (str, optional) – A title for the figure can be provided. (Default:
""
)save (str, optional) – If provided, the figure is saved using
plt.savefig(save)
. A .png or .pdf suffix is recommended. (Default:None
)
- Returns:
out – The preview
Figure
.- Return type:
Figure
- property spectrum
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.
- class martini.martini.Martini(source=None, datacube=None, beam=None, noise=None, sph_kernel=None, spectral_model=None, quiet=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 callinginsert_source_in_cube()
, (optionally)add_noise()
, (optionally)convolve_beam()
andwrite_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.
See also
SPHSource
,DataCube
,martini.beams
,martini.noise
,martini.sph_kernels
,martini.spectral_models
,GlobalProfile
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, velocity_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, channels="velocity") M.write_fits(cubefile, channels="velocity") print(f"Wrote demo fits output to {cubefile}, and beam image to {beamfile}.") try: M.write_hdf5(hdf5file, channels="velocity") except ModuleNotFoundError: print("h5py package not present, skipping hdf5 output demo.") else: print(f"Wrote demo hdf5 output to {hdf5file}.")
- insert_source_in_cube(skip_validation=False, progressbar=None, ncpu=1)[source]
Populates 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 aRuntimeError
if so. This validation can be skipped (at the cost of accuracy!) by setting this parameterTrue
. (Default:False
)progressbar (bool, optional) – A progress bar is shown by default. Progress bars work, with perhaps some visual glitches, in parallel. If
Martini
was initialised withquiet
set toTrue
, progress bars are switched off unless explicitly turned on. (Default:None
)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 asmultiprocessing
). (Default:1
)
- preview(max_points=5000, fig=1, lim=None, vlim=None, point_scaling='auto', title='', save=None)
Produce a figure showing the source particle coordinates and velocities, and the data cube region.
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. (Default:
1000
)fig (int, optional) – Number of the figure in matplotlib, it will be created as
plt.figure(fig)
. (Default:1
)lim (Quantity, optional) –
Quantity
with dimensions of length. The coordinate axes extend from-lim
tolim
. If unspecified, the maximum absolute coordinate 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, the string"datacube"
can be passed. In this case the axis limits are set by the extent of the data cube. (Default:None
)vlim (Quantity, optional) –
Quantity
with dimensions of speed. The velocity axes and colour bar extend from-vlim
tovlim
. 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. (Default:None
)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. (Default:"auto"
)title (str, optional) – A title for the figure can be provided. (Default:
""
)save (str, optional) – If provided, the figure is saved using
plt.savefig(save)
. A .png or .pdf suffix is recommended. (Default:None
)
- Returns:
out – The preview
Figure
.- Return type:
Figure
- write_beam_fits(filename, channels='frequency', overwrite=True)[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.channels (str, optional) – Type of units (
"velocity"
or"frequency"
) used along the spectral axis in output file. (Default:"frequency"
)overwrite (bool, optional) – Whether to allow overwriting existing files. (Default:
True
)
- Raises:
ValueError – If
Martini
was initialized without abeam
.
- write_fits(filename, channels='frequency', overwrite=True)[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.channels (str, optional) – Type of units (
"velocity"
or"frequency"
) used along the spectral axis in output file. (Default:"frequency"
)overwrite (bool, optional) – Whether to allow overwriting existing files. (Default:
True
)
- write_hdf5(filename, channels='frequency', overwrite=True, memmap=False, compact=False)[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.channels (str, optional) – Type of units (
"velocity"
or"frequency"
) used along the spectral axis in output file. (Default:"frequency"
)overwrite (bool, optional) – Whether to allow overwriting existing files. (Default:
True
)memmap (bool, optional) – If
True
, create a file-like object in memory and return it instead of writing file to disk. (Default:False
)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. (Default:False
)