Repository: kyleaoman/martini (language: python)

Issue:
Add a global profile mode
Sometimes you just want a spectrum. Simplified output formats for this case, and simplified/optimised source insertion. Need just a single "pixel" and to calculate the convolution of the beam and sph kernel for each particle. This could be via a new "dataspectrum" mirroring datacube module and new "insert_source_in_spectrum" function mirroring insert_source_in_cube.


Discussion / hints:

Fix patch:
diff --git a/docs/source/global_profile/index.rst b/docs/source/global_profile/index.rst
new file mode 100644
index 0000000..51d9a88
--- /dev/null
+++ b/docs/source/global_profile/index.rst
@@ -0,0 +1,97 @@
+Global profile mode
+===================
+
+Sometimes all you want is a spatially-integrated spectrum of a source. The :class:`~martini.martini.GlobalProfile` class offers a simplified setup of MARTINI 
+
+Assumptions and limitations
+---------------------------
+
+In order to offer a simpler and faster way to produce a spectrum, the :class:`~martini.martini.GlobalProfile` class makes some assumptions:
+
+ - No spatial aperture is assumed. 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 :class:`~martini.martini.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, 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. For example, if the spectrum within a spatial mask defined by a signal-to-noise or other cut is desired, or if spatially-dependent effects like primary beam attenuation are relevant, then the :class:`~martini.martini.GlobalProfile` class should not be used. The :class:`~martini.martini.GlobalProfile` class is mainly intended to efficiently provide a "quick look" at the spectrum, or a reference "ideal" spectrum.
+
+Usage
+-----
+
+The :doc:`source </sources/index>` and :doc:`spectral model </spectral_models/index>` modules should be set up as for a full MARTINI mock observation. The :doc:`beam </beams/index>`, :doc:`noise </noise/index>` and :doc:`sph_kernel </sph_kernels/index>` modules are not relevant. The :doc:`datacube </datacube/index>` module is not used, but a subset of its configuration options are instead given directly to the :class:`~martini.martini.GlobalProfile`. Schematically, an example initialisation looks like:
+
+.. code-block:: python
+
+    from martini import GlobalProfile
+    from martini.sources import SPHSource
+    from martini.spectral_models import GaussianSpectrum
+		
+    source = SPHSource(...)
+    spectral_model = GaussianSpectrum(...)
+
+    gp = GlobalProfile(
+        source=source,
+	spectral_model=spectral_model,
+	n_channels=64,
+	channel_width=10 * U.km * U.s**-1,
+	velocity_centre=source.vsys,
+	channels="velocity",
+    )
+
+The arguments to the other modules are omitted here (replaced with ``...``), check the documentation pages of each module for details. Here the spectrum will be centred on the source systemic velocity, but an explicit frequency or Doppler velocity value could be given instead. The ``channels`` argument determines whether the resulting spectrum will have channel edges in velocity or frequency units. The units (frequency or velocity) of the ``channel_width`` and ``velocity_centre``, and whether ``channels`` is set to ``"frequency"`` or ``"velocity"``, can be mixed in any combination.
+
+Inserting the source
+--------------------
+
+The spectrum can be accessed through the ``spectrum`` attribute:
+
+.. code-block:: python
+
+    gp.spectrum
+
+If it has not yet been calculated, it will be calculated when accessed. The calculation can be explicitly forced with:
+
+.. code-block:: python
+
+    gp.insert_source_in_spectrum()
+
+but this is not usually necessary. In addition to the spectrum itself, the centres and edges of the channels are available as:
+
+.. code-block:: python
+
+    gp.channel_mids
+    gp.channel_edges
+
+respectively. These arrays will have dimensions of frequency or velocity depending on the ``channels`` argument passed to :class:`~martini.martini.GlobalProfile` on initialization.
+
+Parallelization
++++++++++++++++
+
+The core loop in the source insertion function is a loop over pixels. Since parallelization is implemented for this loop, and for a :class:`~martini.martini.GlobalProfile` there is a single pixel, parallelization is not available in this mode.
+
+Quick plot of the spectrum
+--------------------------
+
+As a convenience, a function is provided to make a quick plot showing the spectrum. Whether the systemic velocity of the source (as reported by ``source.vsys``) is shown by a vertical dotted line is controlled by the ``show_vsys`` argument.
+
+.. code-block:: python
+
+    from martini import demo_source, GlobalProfile
+    from martini.spectral_models import GaussianSpectrum
+
+    source = demo_source(N=20000)  # create a simple disc with 20000 particles
+    gp = GlobalProfile(
+        source=source,
+	spectral_model=GaussianSpectrum(sigma=7 * U.km * U.s**-1),
+	n_channels=128,
+	channel_width=2.5 * U.km * U.s**-1,
+	velocity_centre=source.vsys,
+	channels="velocity",
+    )
+    gp.plot_spectrum(
+        show_vsys=True,
+    )
+
+The resulting figure is returned by the function, or can be directly saved to file with the ``save`` argument (e.g. ``save="myplot.png"`` or ``save="myplot.pdf"``). This example looks like:
+
+.. image:: spectrum1.png
+    :width: 350
+    :alt: Example spectrum of a simple rotating disc source.
diff --git a/docs/source/global_profile/spectrum1.png b/docs/source/global_profile/spectrum1.png
new file mode 100644
index 0000000..cc70ec4
Binary files /dev/null and b/docs/source/global_profile/spectrum1.png differ
diff --git a/docs/source/index.rst b/docs/source/index.rst
index d83d667..6701224 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -21,6 +21,7 @@ Overview
    sph_kernels/index
    spectral_models/index
    martini/index
+   global_profile/index
    modules/index
 
 Citing Martini
diff --git a/docs/source/modules/martini.martini.rst b/docs/source/modules/martini.martini.rst
index aa87d6e..c2416ef 100644
--- a/docs/source/modules/martini.martini.rst
+++ b/docs/source/modules/martini.martini.rst
@@ -4,3 +4,4 @@ martini module
 .. automodule:: martini.martini
    :members:
    :show-inheritance:
+   :inherited-members:
diff --git a/martini/__init__.py b/martini/__init__.py
index 4a2e3c1..0f5de86 100644
--- a/martini/__init__.py
+++ b/martini/__init__.py
@@ -1,4 +1,4 @@
-from .martini import Martini
+from .martini import Martini, GlobalProfile
 from .datacube import DataCube
 from ._demo import demo, demo_source
 from .__version__ import __version__
diff --git a/martini/datacube.py b/martini/datacube.py
index fe16063..56b595b 100644
--- a/martini/datacube.py
+++ b/martini/datacube.py
@@ -464,3 +464,42 @@ def __repr__(self):
             Text representation of the DataCube._array contents.
         """
         return self._array.__repr__()
+
+
+class _GlobalProfileDataCube(DataCube):
+    """
+    Helper class that configures a data cube with a single pixel to hold a spectrum.
+
+    Parameters
+    ----------
+    n_channels : int, optional
+        Number of channels along the spectral axis. (Default: 64.)
+
+    channel_width : Quantity, with dimensions of velocity or frequency, optional
+        Step size along the spectral axis. Can be provided as a velocity or a
+        frequency. (Default: 4 km/s.)
+
+    velocity_centre : Quantity, with dimensions of velocity or frequency, optional
+        Velocity (or frequency) of the centre along the spectral axis.
+        (Default: 0 km/s.)
+    """
+
+    def __init__(
+        self,
+        n_channels=64,
+        channel_width=4.0 * U.km * U.s**-1,
+        velocity_centre=0.0 * U.km * U.s**-1,
+    ):
+        super().__init__(
+            n_px_x=1,
+            n_px_y=1,
+            n_channels=n_channels,
+            px_size=1 * U.deg,  # must be >0, ignored for insertion, needed for units
+            channel_width=channel_width,
+            velocity_centre=velocity_centre,
+            ra=0.0 * U.deg,
+            dec=0.0 * U.deg,
+            stokes_axis=False,
+        )
+
+        return
diff --git a/martini/datacube.pyi b/martini/datacube.pyi
index c4fd2ce..ebe578b 100644
--- a/martini/datacube.pyi
+++ b/martini/datacube.pyi
@@ -52,3 +52,11 @@ class DataCube:
     def save_state(self, filename: str, overwrite: bool = ...) -> None: ...
     @classmethod
     def load_state(cls, filename: str) -> "DataCube": ...
+
+class _GlobalProfileDataCube(DataCube):
+    def __init__(
+        self,
+        n_channels: int = ...,
+        channel_width: U.Quantity[U.arcsec] = ...,
+        velocity_centre: U.Quantity[U.km / U.s] = ...,
+    ) -> None: ...
diff --git a/martini/martini.py b/martini/martini.py
index e576401..6efd5ef 100644
--- a/martini/martini.py
+++ b/martini/martini.py
@@ -10,7 +10,8 @@
 from itertools import product
 from .__version__ import __version__ as martini_version
 from warnings import warn
-from martini.datacube import DataCube
+from martini.datacube import DataCube, _GlobalProfileDataCube, HIfreq
+from martini.sph_kernels import DiracDeltaKernel
 
 try:
     gc = subprocess.check_output(
@@ -24,20 +25,13 @@
     martini_version = martini_version + "_commit_" + gc.strip().decode()
 
 
-class Martini:
+class _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.
+    Common methods for the core Martini class and related classes.
 
     Parameters
     ----------
-    source : an instance of a class derived from martini.source._BaseSource
+    source : an instance of a class derived from martini.sources.SPHSource
         A description of the HI emitting object, including position, geometry
         and an interface to the simulation data (SPH particle masses,
         positions, etc.). Sources leveraging the simobj package for reading
@@ -77,107 +71,13 @@ class Martini:
 
     See Also
     --------
-    martini.sources
-    martini.DataCube
+    martini.sources.sph_source.SPHSource
+    martini.datacube.DataCube
     martini.beams
     martini.noise
     martini.sph_kernels
     martini.spectral_models
-
-    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}.")
+    martini.martini.GlobalProfile
     """
 
     def __init__(
@@ -188,6 +88,7 @@ def __init__(
         noise=None,
         sph_kernel=None,
         spectral_model=None,
+        _prune_kwargs=dict(),
         quiet=False,
     ):
         self.quiet = quiet
@@ -196,7 +97,7 @@ def __init__(
         else:
             raise ValueError("A source instance is required.")
         if datacube is not None:
-            self.datacube = datacube
+            self._datacube = datacube
         else:
             raise ValueError("A datacube instance is required.")
         self.beam = beam
@@ -211,97 +112,39 @@ def __init__(
             raise ValueError("A spectral model instance is required.")
 
         if self.beam is not None:
-            self.beam.init_kernel(self.datacube)
-            self.datacube.add_pad(self.beam.needs_pad())
+            self.beam.init_kernel(self._datacube)
+            self._datacube.add_pad(self.beam.needs_pad())
 
         self.source._init_skycoords()
-        self.source._init_pixcoords(self.datacube)  # after datacube is padded
+        self.source._init_pixcoords(self._datacube)  # after datacube is padded
 
-        self.sph_kernel._init_sm_lengths(source=self.source, datacube=self.datacube)
+        self.sph_kernel._init_sm_lengths(source=self.source, datacube=self._datacube)
         self.sph_kernel._init_sm_ranges()
-        self._prune_particles()  # prunes both source, and kernel if applicable
-
-        self.spectral_model.init_spectra(self.source, self.datacube)
-
-        return
-
-    def convolve_beam(self):
-        """
-        Convolve the beam and DataCube.
-        """
-
-        if self.beam is None:
-            warn("Skipping beam convolution, no beam object provided to " "Martini.")
-            return
-
-        minimum_padding = self.beam.needs_pad()
-        if (self.datacube.padx < minimum_padding[0]) or (
-            self.datacube.pady < minimum_padding[1]
-        ):
-            raise ValueError(
-                "datacube padding insufficient for beam convolution (perhaps you loaded a"
-                " datacube state with datacube.load_state that was previously initialized"
-                " by martini with a smaller beam?)"
-            )
-
-        unit = self.datacube._array.unit
-        for spatial_slice in self.datacube.spatial_slices():
-            # use a view [...] to force in-place modification
-            spatial_slice[...] = (
-                fftconvolve(spatial_slice, self.beam.kernel, mode="same") * unit
-            )
-        self.datacube.drop_pad()
-        self.datacube._array = self.datacube._array.to(
-            U.Jy * U.beam**-1,
-            equivalencies=U.beam_angular_area(self.beam.area),
-        )
-        if not self.quiet:
-            print(
-                "Beam convolved.",
-                "  Data cube RMS after beam convolution:"
-                f" {np.std(self.datacube._array):.2e}",
-                f"  Maximum pixel: {self.datacube._array.max():.2e}",
-                "  Median non-zero pixel:"
-                f" {np.median(self.datacube._array[self.datacube._array > 0]):.2e}",
-                sep="\n",
-            )
-        return
-
-    def add_noise(self):
-        """
-        Insert noise into the DataCube.
-        """
+        self._prune_particles(
+            **_prune_kwargs
+        )  # prunes both source, and kernel if applicable
 
-        if self.noise is None:
-            warn("Skipping noise, no noise object provided to Martini.")
-            return
+        self.spectral_model.init_spectra(self.source, self._datacube)
 
-        # this unit conversion means noise can be added before or after source insertion:
-        noise_cube = (
-            self.noise.generate(self.datacube, self.beam)
-            .to(
-                U.Jy * U.arcsec**-2,
-                equivalencies=U.beam_angular_area(self.beam.area),
-            )
-            .to(self.datacube._array.unit, equivalencies=[self.datacube.arcsec2_to_pix])
-        )
-        self.datacube._array = self.datacube._array + noise_cube
-        if not self.quiet:
-            print(
-                "Noise added.",
-                f"  Noise cube RMS: {np.std(noise_cube):.2e} (before beam convolution).",
-                "  Data cube RMS after noise addition (before beam convolution): "
-                f"{np.std(self.datacube._array):.2e}",
-                sep="\n",
-            )
         return
 
-    def _prune_particles(self):
+    def _prune_particles(self, spatial=True, spectral=True, obj_type_str="data cube"):
         """
         Determines which particles cannot contribute to the DataCube and
         removes them to speed up calculation. Assumes the kernel is 0 at
         distances greater than the kernel size (which may differ from the
         SPH smoothing length).
+
+        Parameters
+        ----------
+        spatial : bool
+            If True, prune particles that fall outside the spatial aperture.
+            (Default: True)
+        spectral : bool
+            If True, prune particles that fall outside the spectral bandwidth.
+            (Default: True)
+        obj_type_str : str
+            String describing the object to be pruned for messages. (Default: "data cube")
         """
 
         if not self.quiet:
@@ -310,30 +153,39 @@ def _prune_particles(self):
                 f" mass of {self.source.mHI_g.sum():.2e}."
             )
         spectrum_half_width = (
-            self.spectral_model.half_width(self.source) / self.datacube.channel_width
+            self.spectral_model.half_width(self.source) / self._datacube.channel_width
+        )
+        spatial_reject_conditions = (
+            (
+                (
+                    self.source.pixcoords[:2] + self.sph_kernel.sm_ranges[np.newaxis]
+                    < 0 * U.pix
+                ).any(axis=0),
+                self.source.pixcoords[0] - self.sph_kernel.sm_ranges
+                > (self._datacube.n_px_x + self._datacube.padx * 2) * U.pix,
+                self.source.pixcoords[1] - self.sph_kernel.sm_ranges
+                > (self._datacube.n_px_y + self._datacube.pady * 2) * U.pix,
+            )
+            if spatial
+            else tuple()
         )
-        reject_conditions = (
+        spectral_reject_conditions = (
             (
-                self.source.pixcoords[:2] + self.sph_kernel.sm_ranges[np.newaxis]
-                < 0 * U.pix
-            ).any(axis=0),
-            self.source.pixcoords[0] - self.sph_kernel.sm_ranges
-            > (self.datacube.n_px_x + self.datacube.padx * 2) * U.pix,
-            self.source.pixcoords[1] - self.sph_kernel.sm_ranges
-            > (self.datacube.n_px_y + self.datacube.pady * 2) * U.pix,
-            self.source.pixcoords[2] + 4 * spectrum_half_width * U.pix < 0 * U.pix,
-            self.source.pixcoords[2] - 4 * spectrum_half_width * U.pix
-            > self.datacube.n_channels * U.pix,
+                self.source.pixcoords[2] + 4 * spectrum_half_width * U.pix < 0 * U.pix,
+                self.source.pixcoords[2] - 4 * spectrum_half_width * U.pix
+                > self._datacube.n_channels * U.pix,
+            )
+            if spectral
+            else tuple()
         )
         reject_mask = np.zeros(self.source.pixcoords[0].shape)
-        for condition in reject_conditions:
+        for condition in spatial_reject_conditions + spectral_reject_conditions:
             reject_mask = np.logical_or(reject_mask, condition)
         self.source.apply_mask(np.logical_not(reject_mask))
-        # most kernels ignore this line, but required by AdaptiveKernel
         self.sph_kernel._apply_mask(np.logical_not(reject_mask))
         if not self.quiet:
             print(
-                f"Pruned particles that will not contribute to data cube, "
+                f"Pruned particles that will not contribute to {obj_type_str}, "
                 f"{self.source.npart} particles remaining with total HI mass of "
                 f"{self.source.mHI_g.sum():.2e}."
             )
@@ -348,7 +200,7 @@ def _evaluate_pixel_spectrum(self, ranks_and_ij_pxs, progressbar=True):
         in one-pixel pieces. This avoids the need for concurrent access to the datacube
         by parallel processes, which would in the simplest case duplicate a copy of the
         datacube array per parallel process! In realistic use cases the memory overhead
-        from a the equivalent of a second datacube array should be minimal - memory-
+        from the equivalent of a second datacube array should be minimal - memory-
         limited applications should be limited by the memory consumed by particle data,
         which is not duplicated in parallel execution.
 
@@ -384,7 +236,7 @@ def _evaluate_pixel_spectrum(self, ranks_and_ij_pxs, progressbar=True):
             )
             insertion_slice = (
                 np.s_[ij_px[0], ij_px[1], :, 0]
-                if self.datacube.stokes_axis
+                if self._datacube.stokes_axis
                 else np.s_[ij_px[0], ij_px[1], :]
             )
             result.append(
@@ -405,13 +257,16 @@ def _insert_pixel(self, insertion_slice, insertion_data):
         ----------
         insertion_slice : integer, tuple or slice
             Index into the datacube's _array specifying the insertion location.
+
         insertion data : array-like
             1D array containing the spectrum at the location specified by insertion_slice.
         """
-        self.datacube._array[insertion_slice] = insertion_data
+        self._datacube._array[insertion_slice] = insertion_data
         return
 
-    def insert_source_in_cube(self, skip_validation=False, progressbar=None, ncpu=1):
+    def _insert_source_in_cube(
+        self, skip_validation=False, progressbar=None, ncpu=1, quiet=None
+    ):
         """
         Populates the DataCube with flux from the particles in the source.
 
@@ -436,6 +291,9 @@ def insert_source_in_cube(self, skip_validation=False, progressbar=None, ncpu=1)
             one cpu requires the `multiprocess` module (n.b. not the same as
             `multiprocessing`). (Default: 1)
 
+        quiet : bool, optional
+            If True, suppress output to stdout. If specified, takes precedence over
+            quiet parameter of class. (Default: None)
         """
 
         assert self.spectral_model.spectra is not None
@@ -447,8 +305,8 @@ def insert_source_in_cube(self, skip_validation=False, progressbar=None, ncpu=1)
 
         ij_pxs = list(
             product(
-                np.arange(self.datacube._array.shape[0]),
-                np.arange(self.datacube._array.shape[1]),
+                np.arange(self._datacube._array.shape[0]),
+                np.arange(self._datacube._array.shape[1]),
             )
         )
 
@@ -469,90 +327,637 @@ def insert_source_in_cube(self, skip_validation=False, progressbar=None, ncpu=1)
                     for insertion_slice, insertion_data in result:
                         self._insert_pixel(insertion_slice, insertion_data)
 
-        self.datacube._array = self.datacube._array.to(
-            U.Jy / U.arcsec**2, equivalencies=[self.datacube.arcsec2_to_pix]
+        self._datacube._array = self._datacube._array.to(
+            U.Jy / U.arcsec**2, equivalencies=[self._datacube.arcsec2_to_pix]
         )
         pad_mask = (
             np.s_[
-                self.datacube.padx : -self.datacube.padx,
-                self.datacube.pady : -self.datacube.pady,
+                self._datacube.padx : -self._datacube.padx,
+                self._datacube.pady : -self._datacube.pady,
                 ...,
             ]
-            if self.datacube.padx > 0 and self.datacube.pady > 0
+            if self._datacube.padx > 0 and self._datacube.pady > 0
             else np.s_[...]
         )
-        inserted_flux = self.datacube._array[pad_mask].sum() * self.datacube.px_size**2
+        inserted_flux_density = (
+            self._datacube._array[pad_mask].sum() * self._datacube.px_size**2
+        ).to(U.Jy)
         inserted_mass = (
             2.36e5
             * U.Msun
             * self.source.distance.to_value(U.Mpc) ** 2
-            * inserted_flux.to_value(U.Jy)
-            * self.datacube.channel_width.to_value(U.km / U.s)
+            * inserted_flux_density.to_value(U.Jy)
+            * self._datacube.channel_width.to_value(U.km / U.s)
         )
-        if not self.quiet:
+        if (quiet is None and not self.quiet) or (quiet is not None and not quiet):
             print(
                 "Source inserted.",
-                f"  Flux in cube: {inserted_flux:.2e}",
-                f"  Mass in cube (assuming distance {self.source.distance:.2f}):"
+                f"  Flux density in cube: {inserted_flux_density:.2e}",
+                f"  Mass in cube (assuming distance {self.source.distance:.2f} and a"
+                f" spatially resolved source):"
                 f" {inserted_mass:.2e}",
                 f"    [{inserted_mass / self.source.input_mass * 100:.0f}%"
                 f" of initial source mass]",
-                f"  Maximum pixel: {self.datacube._array.max():.2e}",
+                f"  Maximum pixel: {self._datacube._array.max():.2e}",
                 "  Median non-zero pixel:"
-                f" {np.median(self.datacube._array[self.datacube._array > 0]):.2e}",
+                f" {np.median(self._datacube._array[self._datacube._array > 0]):.2e}",
                 sep="\n",
             )
         return
 
-    def write_fits(
+    def reset(self):
+        """
+        Re-initializes the DataCube with zero-values.
+        """
+        init_kwargs = dict(
+            n_px_x=self._datacube.n_px_x,
+            n_px_y=self._datacube.n_px_y,
+            n_channels=self._datacube.n_channels,
+            px_size=self._datacube.px_size,
+            channel_width=self._datacube.channel_width,
+            velocity_centre=self._datacube.velocity_centre,
+            ra=self._datacube.ra,
+            dec=self._datacube.dec,
+            stokes_axis=self._datacube.stokes_axis,
+        )
+        self._datacube = DataCube(**init_kwargs)
+        if self.beam is not None:
+            self._datacube.add_pad(self.beam.needs_pad())
+        return
+
+    def preview(
         self,
-        filename,
-        channels="frequency",
-        overwrite=True,
+        max_points=5000,
+        fig=1,
+        lim=None,
+        vlim=None,
+        point_scaling="auto",
+        title="",
+        save=None,
     ):
         """
-        Output the DataCube to a FITS-format file.
+        Produce a figure showing the source particle coordinates and velocities, and the
+        datacube 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 datacube extent in position and velocity at the
+        distance and systemic velocity of the source.
 
         Parameters
         ----------
-        filename : string
-            Name of the file to write. '.fits' will be appended if not already
-            present.
+        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)
 
-        channels : {'frequency', 'velocity'}, optional
-            Type of units used along the spectral axis in output file.
-            (Default: 'frequency'.)
+        fig : int, optional
+            Number of the figure in matplotlib, it will be created as `plt.figure(fig)`.
+            (Default: 1)
 
-        overwrite: bool, optional
-            Whether to allow overwriting existing files. (Default: True.)
-        """
+        lim : Quantity with dimensions of length, optional
+            The coordinate axes extend from -lim to lim. 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)
 
-        self.datacube.drop_pad()
-        if channels == "frequency":
-            self.datacube.freq_channels()
-        elif channels == "velocity":
-            self.datacube.velocity_channels()
-        else:
-            raise ValueError(
-                "Martini.write_fits: Unknown 'channels' value "
-                "(use 'frequency' or 'velocity')."
-            )
+        vlim : Quantity with dimensions of speed, optional
+            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, the string "datacube" can be passed. In this
+            case the axis limits are set by the extent of the data cube. (Default: None)
 
-        filename = filename if filename[-5:] == ".fits" else filename + ".fits"
+        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")
 
-        wcs_header = self.datacube.wcs.to_header()
-        wcs_header.rename_keyword("WCSAXES", "NAXIS")
+        title : str, optional
+            A title for the figure can be provided. (Default: "")
 
-        header = fits.Header()
-        header.append(("SIMPLE", "T"))
-        header.append(("BITPIX", 16))
-        header.append(("NAXIS", wcs_header["NAXIS"]))
-        header.append(("NAXIS1", self.datacube.n_px_x))
-        header.append(("NAXIS2", self.datacube.n_px_y))
-        header.append(("NAXIS3", self.datacube.n_channels))
-        if self.datacube.stokes_axis:
-            header.append(("NAXIS4", 1))
-        header.append(("EXTEND", "T"))
+        save : str, optional
+            If provided, the figure is saved using `plt.savefig(save)`. A `.png` or `.pdf`
+            suffix is recommended. (Default: None)
+
+        Returns
+        -------
+        out : matplotlib.figure instance
+            The preview figure.
+        """
+        import matplotlib.pyplot as plt
+
+        y_cube = (
+            (self._datacube.ra - self.source.ra)
+            / np.cos(self.source.dec)
+            * self.source.distance
+        ).to(U.kpc, U.dimensionless_angles())
+        z_cube = ((self._datacube.dec - self.source.dec) * self.source.distance).to(
+            U.kpc, U.dimensionless_angles()
+        )
+        v_cube = (self._datacube.velocity_centre - self.source.vsys).to(U.km / U.s)
+        dy_cube = 0.5 * (
+            self._datacube.px_size * self._datacube.n_px_x * self.source.distance
+        ).to(
+            U.kpc, U.dimensionless_angles()
+        )  # half-length
+        dz_cube = 0.5 * (
+            self._datacube.px_size * self._datacube.n_px_y * self.source.distance
+        ).to(
+            U.kpc, U.dimensionless_angles()
+        )  # half-length
+        dv_cube = 0.5 * (self._datacube.channel_width * self._datacube.n_channels).to(
+            U.km / U.s
+        )
+        # should issue some warnings if dRA or dDec > 5deg, or if dec within 5 deg of pole
+        clip_lim = lim == "datacube"
+        if clip_lim:
+            lim = max(
+                np.abs(
+                    U.Quantity(
+                        [
+                            y_cube - dy_cube,
+                            y_cube + dy_cube,
+                            z_cube - dz_cube,
+                            z_cube + dz_cube,
+                        ]
+                    )
+                )
+            )
+        clip_vlim = vlim == "datacube"
+        if clip_vlim:
+            clip_vlim = True
+            vlim = max(np.abs(U.Quantity([v_cube - dv_cube, v_cube + dv_cube])))
+        else:
+            clip_vlim = False
+        # pass through arguments, except save (which we will do later if desired)
+        fig = self.source.preview(
+            max_points=max_points,
+            fig=fig,
+            lim=lim,
+            vlim=vlim,
+            point_scaling=point_scaling,
+            title=title,
+            save=None,
+        )
+        sp1, cb, sp2, sp3 = fig.get_axes()
+        sp1.add_patch(
+            plt.Rectangle(
+                (
+                    (y_cube - dy_cube).to_value(U.kpc),
+                    (z_cube - dz_cube).to_value(U.kpc),
+                ),
+                2 * dy_cube.to_value(U.kpc),
+                2 * dz_cube.to_value(U.kpc),
+                linewidth=5,
+                edgecolor="red",
+                facecolor="None",
+            )
+        )
+        sp2.add_patch(
+            plt.Rectangle(
+                (
+                    (y_cube - dy_cube).to_value(U.kpc),
+                    (v_cube - dv_cube).to_value(U.km / U.s),
+                ),
+                2 * dy_cube.to_value(U.kpc),
+                2 * dv_cube.to_value(U.km / U.s),
+                linewidth=5,
+                edgecolor="red",
+                facecolor="None",
+            )
+        )
+        sp3.add_patch(
+            plt.Rectangle(
+                (
+                    (z_cube - dz_cube).to_value(U.kpc),
+                    (v_cube - dv_cube).to_value(U.km / U.s),
+                ),
+                2 * dz_cube.to_value(U.kpc),
+                2 * dv_cube.to_value(U.km / U.s),
+                linewidth=5,
+                edgecolor="red",
+                facecolor="None",
+            )
+        )
+        if clip_lim:
+            sp1.set_xlim(
+                (y_cube + dy_cube).to_value(U.kpc), (y_cube - dy_cube).to_value(U.kpc)
+            )
+            sp1.set_ylim(
+                (z_cube - dz_cube).to_value(U.kpc), (z_cube + dz_cube).to_value(U.kpc)
+            )
+            sp2.set_xlim(
+                (y_cube + dy_cube).to_value(U.kpc), (y_cube - dy_cube).to_value(U.kpc)
+            )
+            sp3.set_xlim(
+                (z_cube - dz_cube).to_value(U.kpc), (z_cube + dz_cube).to_value(U.kpc)
+            )
+        else:
+            if lim is None:
+                lim = sp1.get_xlim()[1] * U.kpc
+            sp1.set_xlim(
+                (
+                    max(1.1 * (y_cube + dy_cube), lim).to_value(U.kpc),
+                    min(1.1 * (y_cube - dy_cube), -lim).to_value(U.kpc),
+                )
+            )
+            sp2.set_xlim(
+                (
+                    max(1.1 * (y_cube + dy_cube), lim).to_value(U.kpc),
+                    min(1.1 * (y_cube - dy_cube), -lim).to_value(U.kpc),
+                )
+            )
+            sp1.set_ylim(
+                (
+                    min(1.1 * (z_cube - dz_cube), -lim).to_value(U.kpc),
+                    max(1.1 * (z_cube + dz_cube), lim).to_value(U.kpc),
+                )
+            )
+            sp3.set_xlim(
+                (
+                    min(1.1 * (z_cube - dz_cube), -lim).to_value(U.kpc),
+                    max(1.1 * (z_cube + dz_cube), lim).to_value(U.kpc),
+                )
+            )
+        if clip_vlim:
+            sp2.set_ylim(
+                (v_cube - dv_cube).to_value(U.km / U.s),
+                (v_cube + dv_cube).to_value(U.km / U.s),
+            )
+            sp3.set_ylim(
+                (v_cube - dv_cube).to_value(U.km / U.s),
+                (v_cube + dv_cube).to_value(U.km / U.s),
+            )
+        else:
+            if vlim is None:
+                vlim = sp2.get_ylim()[1] * U.km / U.s
+            sp2.set_ylim(
+                (
+                    min(1.1 * (v_cube - dv_cube), -vlim).to_value(U.km / U.s),
+                    max(1.1 * (v_cube + dv_cube), vlim).to_value(U.km / U.s),
+                )
+            )
+            sp3.set_ylim(
+                (
+                    min(1.1 * (v_cube - dv_cube), -vlim).to_value(U.km / U.s),
+                    max(1.1 * (v_cube + dv_cube), vlim).to_value(U.km / U.s),
+                )
+            )
+
+        if save is not None:
+            plt.savefig(save)
+        return fig
+
+
+class Martini(_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 : an instance of a class derived from martini.sources.sph_source.SPHSource
+        A description of the HI emitting object, including position, geometry
+        and an interface to the simulation data (SPH particle masses,
+        positions, etc.). Sources leveraging the simobj package for reading
+        simulation data (github.com/kyleaoman/simobj) and a few test sources
+        (e.g. single particle) are provided, creation of customized sources,
+        for instance to leverage other interfaces to simulation data, is
+        straightforward. See sub-module documentation.
+
+    datacube : martini.DataCube instance
+        A description of the datacube to create, including pixels, channels,
+        sky position. See sub-module documentation.
+
+    beam : an instance of a class derived from beams._BaseBeam, optional
+        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 : an instance of a class derived from noise._BaseNoise, optional
+        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 : an instance of a class derived from sph_kernels._BaseSPHKernel
+        A description of the SPH smoothing kernel. Check simulation
+        documentation for the kernel used in a particular simulation, and
+        SPH kernel submodule documentation for guidance.
+
+    spectral_model : an instance of a class derived from \
+    spectral_models._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
+        If True, suppress output to stdout. (Default: False)
+
+    See Also
+    --------
+    martini.sources.sph_source.SPHSource
+    martini.datacube.DataCube
+    martini.beams
+    martini.noise
+    martini.sph_kernels
+    martini.spectral_models
+    martini.martini.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}.")
+    """
+
+    def __init__(
+        self,
+        source=None,
+        datacube=None,
+        beam=None,
+        noise=None,
+        sph_kernel=None,
+        spectral_model=None,
+        quiet=False,
+    ):
+        super().__init__(
+            source=source,
+            datacube=datacube,
+            beam=beam,
+            noise=noise,
+            sph_kernel=sph_kernel,
+            spectral_model=spectral_model,
+            quiet=quiet,
+        )
+
+        return
+
+    @property
+    def datacube(self):
+        """
+        The `martini.datacube.DataCube` object for this mock observation.
+        """
+        return self._datacube
+
+    def insert_source_in_cube(self, skip_validation=False, progressbar=None, ncpu=1):
+        """
+        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 a
+            RuntimeError if so. This validation can be skipped (at the cost
+            of accuracy!) by setting this parameter True. (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 with
+            `quiet` set to `True`, 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 as
+            `multiprocessing`). (Default: 1)
+
+        """
+
+        super()._insert_source_in_cube(
+            skip_validation=skip_validation, progressbar=progressbar, ncpu=ncpu
+        )
+
+        return
+
+    def convolve_beam(self):
+        """
+        Convolve the beam and DataCube.
+        """
+
+        if self.beam is None:
+            warn("Skipping beam convolution, no beam object provided to " "Martini.")
+            return
+
+        minimum_padding = self.beam.needs_pad()
+        if (self._datacube.padx < minimum_padding[0]) or (
+            self._datacube.pady < minimum_padding[1]
+        ):
+            raise ValueError(
+                "datacube padding insufficient for beam convolution (perhaps you loaded a"
+                " datacube state with datacube.load_state that was previously initialized"
+                " by martini with a smaller beam?)"
+            )
+
+        unit = self._datacube._array.unit
+        for spatial_slice in self._datacube.spatial_slices():
+            # use a view [...] to force in-place modification
+            spatial_slice[...] = (
+                fftconvolve(spatial_slice, self.beam.kernel, mode="same") * unit
+            )
+        self._datacube.drop_pad()
+        self._datacube._array = self._datacube._array.to(
+            U.Jy * U.beam**-1,
+            equivalencies=U.beam_angular_area(self.beam.area),
+        )
+        if not self.quiet:
+            print(
+                "Beam convolved.",
+                "  Data cube RMS after beam convolution:"
+                f" {np.std(self._datacube._array):.2e}",
+                f"  Maximum pixel: {self._datacube._array.max():.2e}",
+                "  Median non-zero pixel:"
+                f" {np.median(self._datacube._array[self._datacube._array > 0]):.2e}",
+                sep="\n",
+            )
+        return
+
+    def add_noise(self):
+        """
+        Insert noise into the DataCube.
+        """
+
+        if self.noise is None:
+            warn("Skipping noise, no noise object provided to Martini.")
+            return
+
+        # this unit conversion means noise can be added before or after source insertion:
+        noise_cube = (
+            self.noise.generate(self._datacube, self.beam)
+            .to(
+                U.Jy * U.arcsec**-2,
+                equivalencies=U.beam_angular_area(self.beam.area),
+            )
+            .to(
+                self._datacube._array.unit,
+                equivalencies=[self._datacube.arcsec2_to_pix],
+            )
+        )
+        self._datacube._array = self._datacube._array + noise_cube
+        if not self.quiet:
+            print(
+                "Noise added.",
+                f"  Noise cube RMS: {np.std(noise_cube):.2e} (before beam convolution).",
+                "  Data cube RMS after noise addition (before beam convolution): "
+                f"{np.std(self._datacube._array):.2e}",
+                sep="\n",
+            )
+        return
+
+    def write_fits(
+        self,
+        filename,
+        channels="frequency",
+        overwrite=True,
+    ):
+        """
+        Output the DataCube to a FITS-format file.
+
+        Parameters
+        ----------
+        filename : string
+            Name of the file to write. '.fits' will be appended if not already
+            present.
+
+        channels : {'frequency', 'velocity'}, optional
+            Type of units used along the spectral axis in output file.
+            (Default: 'frequency'.)
+
+        overwrite: bool, optional
+            Whether to allow overwriting existing files. (Default: True.)
+        """
+
+        self._datacube.drop_pad()
+        if channels == "frequency":
+            self._datacube.freq_channels()
+        elif channels == "velocity":
+            self._datacube.velocity_channels()
+        else:
+            raise ValueError(
+                "Martini.write_fits: Unknown 'channels' value "
+                "(use 'frequency' or 'velocity')."
+            )
+
+        filename = filename if filename[-5:] == ".fits" else filename + ".fits"
+
+        wcs_header = self._datacube.wcs.to_header()
+        wcs_header.rename_keyword("WCSAXES", "NAXIS")
+
+        header = fits.Header()
+        header.append(("SIMPLE", "T"))
+        header.append(("BITPIX", 16))
+        header.append(("NAXIS", wcs_header["NAXIS"]))
+        header.append(("NAXIS1", self._datacube.n_px_x))
+        header.append(("NAXIS2", self._datacube.n_px_y))
+        header.append(("NAXIS3", self._datacube.n_channels))
+        if self._datacube.stokes_axis:
+            header.append(("NAXIS4", 1))
+        header.append(("EXTEND", "T"))
         header.append(("CDELT1", wcs_header["CDELT1"]))
         header.append(("CRPIX1", wcs_header["CRPIX1"]))
         header.append(("CRVAL1", wcs_header["CRVAL1"]))
@@ -568,7 +973,7 @@ def write_fits(
         header.append(("CRVAL3", wcs_header["CRVAL3"]))
         header.append(("CTYPE3", wcs_header["CTYPE3"]))
         header.append(("CUNIT3", wcs_header["CUNIT3"]))
-        if self.datacube.stokes_axis:
+        if self._datacube.stokes_axis:
             header.append(("CDELT4", wcs_header["CDELT4"]))
             header.append(("CRPIX4", wcs_header["CRPIX4"]))
             header.append(("CRVAL4", wcs_header["CRVAL4"]))
@@ -579,12 +984,12 @@ def write_fits(
         # header.append(('BLANK', -32768)) #only for integer data
         header.append(("BSCALE", 1.0))
         header.append(("BZERO", 0.0))
-        datacube_array_units = self.datacube._array.unit
+        datacube_array_units = self._datacube._array.unit
         header.append(
-            ("DATAMAX", np.max(self.datacube._array.to_value(datacube_array_units)))
+            ("DATAMAX", np.max(self._datacube._array.to_value(datacube_array_units)))
         )
         header.append(
-            ("DATAMIN", np.min(self.datacube._array.to_value(datacube_array_units)))
+            ("DATAMIN", np.min(self._datacube._array.to_value(datacube_array_units)))
         )
         header.append(("ORIGIN", "astropy v" + astropy_version))
         # long names break fits format, don't let the user set this
@@ -611,12 +1016,12 @@ def write_fits(
 
         # flip axes to write
         hdu = fits.PrimaryHDU(
-            header=header, data=self.datacube._array.to_value(datacube_array_units).T
+            header=header, data=self._datacube._array.to_value(datacube_array_units).T
         )
         hdu.writeto(filename, overwrite=overwrite)
 
         if channels == "frequency":
-            self.datacube.velocity_channels()
+            self._datacube.velocity_channels()
         return
 
     def write_beam_fits(self, filename, channels="frequency", overwrite=True):
@@ -651,9 +1056,9 @@ def write_beam_fits(self, filename, channels="frequency", overwrite=True):
             )
         assert self.beam.kernel is not None
         if channels == "frequency":
-            self.datacube.freq_channels()
+            self._datacube.freq_channels()
         elif channels == "velocity":
-            self.datacube.velocity_channels()
+            self._datacube.velocity_channels()
         else:
             raise ValueError(
                 "Martini.write_beam_fits: Unknown 'channels' "
@@ -662,7 +1067,7 @@ def write_beam_fits(self, filename, channels="frequency", overwrite=True):
 
         filename = filename if filename[-5:] == ".fits" else filename + ".fits"
 
-        wcs_header = self.datacube.wcs.to_header()
+        wcs_header = self._datacube.wcs.to_header()
 
         beam_kernel_units = self.beam.kernel.unit
         header = fits.Header()
@@ -715,7 +1120,7 @@ def write_beam_fits(self, filename, channels="frequency", overwrite=True):
         hdu.writeto(filename, overwrite=True)
 
         if channels == "frequency":
-            self.datacube.velocity_channels()
+            self._datacube.velocity_channels()
         return
 
     def write_hdf5(
@@ -755,9 +1160,9 @@ def write_hdf5(
 
         import h5py
 
-        self.datacube.drop_pad()
+        self._datacube.drop_pad()
         if channels == "frequency":
-            self.datacube.freq_channels()
+            self._datacube.freq_channels()
         elif channels == "velocity":
             pass
         else:
@@ -768,22 +1173,22 @@ def write_hdf5(
 
         filename = filename if filename[-5:] == ".hdf5" else filename + ".hdf5"
 
-        wcs_header = self.datacube.wcs.to_header()
+        wcs_header = self._datacube.wcs.to_header()
 
         mode = "w" if overwrite else "x"
         driver = "core" if memmap else None
         h5_kwargs = {"backing_store": False} if memmap else dict()
         f = h5py.File(filename, mode, driver=driver, **h5_kwargs)
-        datacube_array_units = self.datacube._array.unit
-        s = np.s_[..., 0] if self.datacube.stokes_axis else np.s_[...]
-        f["FluxCube"] = self.datacube._array.to_value(datacube_array_units)[s]
+        datacube_array_units = self._datacube._array.unit
+        s = np.s_[..., 0] if self._datacube.stokes_axis else np.s_[...]
+        f["FluxCube"] = self._datacube._array.to_value(datacube_array_units)[s]
         c = f["FluxCube"]
         origin = 0  # index from 0 like numpy, not from 1
         if not compact:
             xgrid, ygrid, vgrid = np.meshgrid(
-                np.arange(self.datacube._array.shape[0]),
-                np.arange(self.datacube._array.shape[1]),
-                np.arange(self.datacube._array.shape[2]),
+                np.arange(self._datacube._array.shape[0]),
+                np.arange(self._datacube._array.shape[1]),
+                np.arange(self._datacube._array.shape[2]),
             )
             cgrid = (
                 np.vstack(
@@ -794,7 +1199,7 @@ def write_hdf5(
                         np.zeros(vgrid.shape).flatten(),
                     )
                 ).T
-                if self.datacube.stokes_axis
+                if self._datacube.stokes_axis
                 else np.vstack(
                     (
                         xgrid.flatten(),
@@ -803,10 +1208,10 @@ def write_hdf5(
                     )
                 ).T
             )
-            wgrid = self.datacube.wcs.all_pix2world(cgrid, origin)
-            ragrid = wgrid[:, 0].reshape(self.datacube._array.shape)[s]
-            decgrid = wgrid[:, 1].reshape(self.datacube._array.shape)[s]
-            chgrid = wgrid[:, 2].reshape(self.datacube._array.shape)[s]
+            wgrid = self._datacube.wcs.all_pix2world(cgrid, origin)
+            ragrid = wgrid[:, 0].reshape(self._datacube._array.shape)[s]
+            decgrid = wgrid[:, 1].reshape(self._datacube._array.shape)[s]
+            chgrid = wgrid[:, 2].reshape(self._datacube._array.shape)[s]
             f["RA"] = ragrid
             f["RA"].attrs["Unit"] = wcs_header["CUNIT1"]
             f["Dec"] = decgrid
@@ -814,7 +1219,7 @@ def write_hdf5(
             f["channel_mids"] = chgrid
             f["channel_mids"].attrs["Unit"] = wcs_header["CUNIT3"]
         c.attrs["AxisOrder"] = "(RA,Dec,Channels)"
-        c.attrs["FluxCubeUnit"] = str(self.datacube._array.unit)
+        c.attrs["FluxCubeUnit"] = str(self._datacube._array.unit)
         c.attrs["deltaRA_in_RAUnit"] = wcs_header["CDELT1"]
         c.attrs["RA0_in_px"] = wcs_header["CRPIX1"] - 1
         c.attrs["RA0_in_RAUnit"] = wcs_header["CRVAL1"]
@@ -869,91 +1274,298 @@ def write_hdf5(
             b.attrs["MartiniVersion"] = martini_version
             b.attrs["AstropyVersion"] = astropy_version
 
-        if channels == "frequency":
-            self.datacube.velocity_channels()
-        if memmap:
-            return f
+        if channels == "frequency":
+            self._datacube.velocity_channels()
+        if memmap:
+            return f
+        else:
+            f.close()
+            return
+
+
+class GlobalProfile(_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 : an instance of a class derived from martini.sources.sph_source.SPHSource
+        A description of the HI emitting object, including position, geometry
+        and an interface to the simulation data (SPH particle masses,
+        positions, etc.). Sources leveraging the simobj package for reading
+        simulation data (github.com/kyleaoman/simobj) and a few test sources
+        (e.g. single particle) are provided, creation of customized sources,
+        for instance to leverage other interfaces to simulation data, is
+        straightforward. See sub-module documentation.
+
+    spectral_model : an instance of a class derived from \
+    spectral_models._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, with dimensions of velocity or frequency, optional
+        Step size along the spectral axis. Can be provided as a velocity or a
+        frequency. (Default: 4 km/s.)
+
+    velocity_centre : Quantity, with dimensions of velocity or frequency, optional
+        Velocity (or frequency) of the centre along the spectral axis.
+        (Default: 0 km/s.)
+
+    channels : {"frequency", "velocity"}, optional
+            Type of units 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
+    --------
+    martini.sources.sph_source.SPHSource
+    martini.spectral_models
+    martini.martini.Martini
+
+    """
+
+    def __init__(
+        self,
+        source=None,
+        spectral_model=None,
+        n_channels=64,
+        channel_width=4 * U.km * U.s**-1,
+        velocity_centre=0 * U.km * U.s**-1,
+        channels="frequency",
+        quiet=False,
+    ):
+        super().__init__(
+            source=source,
+            datacube=_GlobalProfileDataCube(
+                n_channels=n_channels,
+                channel_width=channel_width,
+                velocity_centre=velocity_centre,
+            ),
+            beam=None,
+            noise=None,
+            sph_kernel=DiracDeltaKernel(size_in_fwhm=np.inf),
+            spectral_model=spectral_model,
+            _prune_kwargs=dict(
+                spatial=False,
+                spectral=True,
+                obj_type_str="spectrum",
+            ),
+            quiet=quiet,
+        )
+        if channels == "velocity":
+            self._datacube.velocity_channels()
+        elif channels == "frequency":
+            self._datacube.freq_channels()
         else:
-            f.close()
-            return
+            raise ValueError(
+                'GlobalProfile: `channels` must be "velocity" or "frequency".'
+            )
+        self.channels = channels
+        self.source.pixcoords[:2] = 0
 
-    def reset(self):
+        return
+
+    def insert_source_in_spectrum(self):
         """
-        Re-initializes the DataCube with zero-values.
+        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.
         """
-        init_kwargs = dict(
-            n_px_x=self.datacube.n_px_x,
-            n_px_y=self.datacube.n_px_y,
-            n_channels=self.datacube.n_channels,
-            px_size=self.datacube.px_size,
-            channel_width=self.datacube.channel_width,
-            velocity_centre=self.datacube.velocity_centre,
-            ra=self.datacube.ra,
-            dec=self.datacube.dec,
-            stokes_axis=self.datacube.stokes_axis,
+        # skip_validation=True: all particles can contribute their kernel to the pixel;
+        # ncpu=1 since we have 1 pixel and source insertion is parallel over pixels;
+        # no progressbar since there's only 1 pixel of progress;
+        # quiet=True because messages assume a resolved source, replace with new ones
+        super()._insert_source_in_cube(
+            skip_validation=True, progressbar=False, ncpu=1, quiet=True
         )
-        self.datacube = DataCube(**init_kwargs)
-        if self.beam is not None:
-            self.datacube.add_pad(self.beam.needs_pad())
+        # The datacube in Jy/arcsec^2 is a bit misleading because the source is
+        # (presumably) completely unresolved so extrapolating its surface brightness
+        # across the entire pixel is incorrect. Correctly integrate out spatial
+        # information and convert to Jy:
+        self._spectrum = (
+            (self._datacube._array.squeeze()).to(
+                U.Jy / U.pix**2, equivalencies=[self._datacube.arcsec2_to_pix]
+            )
+            * U.pix**2
+        ).to(U.Jy)
+        if not self.quiet:
+            # Need a slightly different calculation for a completely unresolved source.
+            inserted_flux_density = self.spectrum.sum()
+            inserted_mass = (
+                2.36e5
+                * U.Msun
+                * self.source.distance.to_value(U.Mpc) ** 2
+                * inserted_flux_density.to_value(U.Jy)
+                * self._datacube.channel_width.to_value(U.km / U.s)
+            )
+            print(
+                "Source inserted.",
+                f"  Flux density in spectrum: {inserted_flux_density:.2e}",
+                f"  Mass in spectrum (assuming distance {self.source.distance:.2f}):"
+                f" {inserted_mass:.2e}",
+                f"    [{inserted_mass / self.source.input_mass * 100:.0f}%"
+                f" of initial source mass]",
+                sep="\n",
+            )
+
+    @property
+    def spectrum(self):
+        """
+        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.
+        """
+        if not hasattr(self, "_spectrum"):
+            self.insert_source_in_spectrum()
+        return self._spectrum
+
+    @property
+    def channel_edges(self):
+        """
+        The edges of the channels for the spectrum.
+        """
+        return self._datacube.channel_edges
+
+    @property
+    def channel_mids(self):
+        """
+        The centres of the channels for the spectrum.
+        """
+        return self._datacube.channel_mids
+
+    @property
+    def channel_width(self):
+        """
+        The width of the channels for the spectrum.
+        """
+        return self._datacube.channel_width
+
+    def reset(self):
+        super().reset()
+        if hasattr(self, "_spectrum"):
+            del self._spectrum
         return
 
-    def preview(
+    def plot_spectrum(
         self,
-        max_points=5000,
         fig=1,
-        lim=None,
-        vlim=None,
-        point_scaling="auto",
         title="",
+        show_vsys=True,
         save=None,
     ):
         """
-        Produce a figure showing the source particle coordinates and velocities, and the
-        datacube 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 datacube extent in position and velocity at the
-        distance and systemic velocity of the source.
+        Produce a figure showing the spectrum.
 
         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 with dimensions of length, optional
-            The coordinate axes extend from -lim to lim. 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 with dimensions of speed, optional
-            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, the string "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: "")
 
+        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)
@@ -961,167 +1573,36 @@ def preview(
         Returns
         -------
         out : matplotlib.figure instance
-            The preview figure.
+            The spectrum figure.
         """
         import matplotlib.pyplot as plt
 
-        y_cube = (
-            (self.datacube.ra - self.source.ra)
-            / np.cos(self.source.dec)
-            * self.source.distance
-        ).to(U.kpc, U.dimensionless_angles())
-        z_cube = ((self.datacube.dec - self.source.dec) * self.source.distance).to(
-            U.kpc, U.dimensionless_angles()
-        )
-        v_cube = (self.datacube.velocity_centre - self.source.vsys).to(U.km / U.s)
-        dy_cube = 0.5 * (
-            self.datacube.px_size * self.datacube.n_px_x * self.source.distance
-        ).to(
-            U.kpc, U.dimensionless_angles()
-        )  # half-length
-        dz_cube = 0.5 * (
-            self.datacube.px_size * self.datacube.n_px_y * self.source.distance
-        ).to(
-            U.kpc, U.dimensionless_angles()
-        )  # half-length
-        dv_cube = 0.5 * (self.datacube.channel_width * self.datacube.n_channels).to(
-            U.km / U.s
-        )
-        # should issue some warnings if dRA or dDec > 5deg, or if dec within 5 deg of pole
-        clip_lim = lim == "datacube"
-        if clip_lim:
-            lim = max(
-                np.abs(
-                    U.Quantity(
-                        [
-                            y_cube - dy_cube,
-                            y_cube + dy_cube,
-                            z_cube - dz_cube,
-                            z_cube + dz_cube,
-                        ]
-                    )
-                )
-            )
-        clip_vlim = vlim == "datacube"
-        if clip_vlim:
-            clip_vlim = True
-            vlim = max(np.abs(U.Quantity([v_cube - dv_cube, v_cube + dv_cube])))
-        else:
-            clip_vlim = False
-        # pass through arguments, except save (which we will do later if desired)
-        fig = self.source.preview(
-            max_points=max_points,
-            fig=fig,
-            lim=lim,
-            vlim=vlim,
-            point_scaling=point_scaling,
-            title=title,
-            save=None,
-        )
-        sp1, cb, sp2, sp3 = fig.get_axes()
-        sp1.add_patch(
-            plt.Rectangle(
-                (
-                    (y_cube - dy_cube).to_value(U.kpc),
-                    (z_cube - dz_cube).to_value(U.kpc),
-                ),
-                2 * dy_cube.to_value(U.kpc),
-                2 * dz_cube.to_value(U.kpc),
-                linewidth=5,
-                edgecolor="red",
-                facecolor="None",
-            )
-        )
-        sp2.add_patch(
-            plt.Rectangle(
-                (
-                    (y_cube - dy_cube).to_value(U.kpc),
-                    (v_cube - dv_cube).to_value(U.km / U.s),
-                ),
-                2 * dy_cube.to_value(U.kpc),
-                2 * dv_cube.to_value(U.km / U.s),
-                linewidth=5,
-                edgecolor="red",
-                facecolor="None",
-            )
-        )
-        sp3.add_patch(
-            plt.Rectangle(
-                (
-                    (z_cube - dz_cube).to_value(U.kpc),
-                    (v_cube - dv_cube).to_value(U.km / U.s),
-                ),
-                2 * dz_cube.to_value(U.kpc),
-                2 * dv_cube.to_value(U.km / U.s),
-                linewidth=5,
-                edgecolor="red",
-                facecolor="None",
-            )
+        fig = plt.figure(fig, figsize=(4, 3))
+        fig.clf()
+        fig.suptitle(title)
+
+        xunit = dict(velocity=U.km * U.s**-1, frequency=U.GHz)[self.channels]
+
+        ax = fig.add_subplot(1, 1, 1)
+        ax.plot(
+            self.channel_mids.to_value(xunit),
+            self.spectrum.to_value(U.Jy),
+            ls="solid",
+            color="black",
+            lw=2,
         )
-        if clip_lim:
-            sp1.set_xlim(
-                (y_cube + dy_cube).to_value(U.kpc), (y_cube - dy_cube).to_value(U.kpc)
-            )
-            sp1.set_ylim(
-                (z_cube - dz_cube).to_value(U.kpc), (z_cube + dz_cube).to_value(U.kpc)
-            )
-            sp2.set_xlim(
-                (y_cube + dy_cube).to_value(U.kpc), (y_cube - dy_cube).to_value(U.kpc)
-            )
-            sp3.set_xlim(
-                (z_cube - dz_cube).to_value(U.kpc), (z_cube + dz_cube).to_value(U.kpc)
-            )
-        else:
-            if lim is None:
-                lim = sp1.get_xlim()[1] * U.kpc
-            sp1.set_xlim(
-                (
-                    max(1.1 * (y_cube + dy_cube), lim).to_value(U.kpc),
-                    min(1.1 * (y_cube - dy_cube), -lim).to_value(U.kpc),
-                )
-            )
-            sp2.set_xlim(
-                (
-                    max(1.1 * (y_cube + dy_cube), lim).to_value(U.kpc),
-                    min(1.1 * (y_cube - dy_cube), -lim).to_value(U.kpc),
-                )
-            )
-            sp1.set_ylim(
-                (
-                    min(1.1 * (z_cube - dz_cube), -lim).to_value(U.kpc),
-                    max(1.1 * (z_cube + dz_cube), lim).to_value(U.kpc),
-                )
-            )
-            sp3.set_xlim(
-                (
-                    min(1.1 * (z_cube - dz_cube), -lim).to_value(U.kpc),
-                    max(1.1 * (z_cube + dz_cube), lim).to_value(U.kpc),
-                )
-            )
-        if clip_vlim:
-            sp2.set_ylim(
-                (v_cube - dv_cube).to_value(U.km / U.s),
-                (v_cube + dv_cube).to_value(U.km / U.s),
-            )
-            sp3.set_ylim(
-                (v_cube - dv_cube).to_value(U.km / U.s),
-                (v_cube + dv_cube).to_value(U.km / U.s),
-            )
-        else:
-            if vlim is None:
-                vlim = sp2.get_ylim()[1] * U.km / U.s
-            sp2.set_ylim(
-                (
-                    min(1.1 * (v_cube - dv_cube), -vlim).to_value(U.km / U.s),
-                    max(1.1 * (v_cube + dv_cube), vlim).to_value(U.km / U.s),
-                )
-            )
-            sp3.set_ylim(
-                (
-                    min(1.1 * (v_cube - dv_cube), -vlim).to_value(U.km / U.s),
-                    max(1.1 * (v_cube + dv_cube), vlim).to_value(U.km / U.s),
-                )
+        if show_vsys:
+            ax.axvline(
+                self.source.vsys.to_value(xunit, equivalencies=U.doppler_radio(HIfreq)),
+                linestyle="dotted",
+                lw=1.5,
+                color="black",
             )
+        ax.set_ylabel(r"Flux density $[\mathrm{Jy}]$")
+        if self.channels == "velocity":
+            ax.set_xlabel(r"Velocity $[\mathrm{km}\,\mathrm{s}^{-1}]$")
+        elif self.channels == "frequency":
+            ax.set_xlabel(r"Frequency $[\mathrm{GHz}]$")
 
         if save is not None:
             plt.savefig(save)
diff --git a/martini/martini.pyi b/martini/martini.pyi
index d64e53b..8561ce4 100644
--- a/martini/martini.pyi
+++ b/martini/martini.pyi
@@ -1,3 +1,4 @@
+from numpy import ndarray
 import typing as T
 from martini.beams import _BaseBeam
 from martini.datacube import DataCube as DataCube
@@ -10,9 +11,9 @@ import astropy.units as U
 
 gc: bytes
 
-class Martini:
+class _BaseMartini:
     source: SPHSource
-    datacube: DataCube
+    _datacube: DataCube
     beam: _BaseBeam
     noise: _BaseNoise
     sph_kernel: _BaseSPHKernel
@@ -27,16 +28,60 @@ class Martini:
         noise: T.Optional[_BaseNoise] = ...,
         sph_kernel: T.Optional[_BaseSPHKernel] = ...,
         spectral_model: T.Optional[_BaseSpectrum] = ...,
+        _prune_kwargs: T.Dict[str, T.Union[bool, str]] = ...,
         quiet: T.Optional[bool] = ...,
     ) -> None: ...
-    def convolve_beam(self) -> None: ...
-    def add_noise(self) -> None: ...
+    def _prune_particles(
+        self, spatial: bool = ..., spectral: bool = ..., obj_type_str: str = ...
+    ) -> None: ...
+    def _evaluate_pixel_spectrum(
+        self,
+        ranks_and_ij_pxs: T.Tuple[int, T.List[T.Tuple[int, int]]],
+        progressbar: bool = ...,
+    ) -> T.Tuple[slice, U.Quantity[U.Jy / U.arcsec**2]]: ...
+    def _insert_pixel(
+        self, insertion_slice: T.Union[int, T.Tuple, slice], insertion_data: ndarray
+    ) -> None: ...
+    def _insert_source_in_cube(
+        self,
+        skip_validation: bool = ...,
+        progressbar: T.Optional[bool] = ...,
+        ncpu: int = ...,
+        quiet: T.Optional[bool] = ...,
+    ) -> None: ...
+    def reset(self) -> None: ...
+    def preview(
+        self,
+        max_points: int = ...,
+        fig: int = ...,
+        lim: T.Optional[T.Union[str, U.Quantity[U.kpc]]] = ...,
+        vlim: T.Optional[T.Union[str, U.Quantity[U.km / U.s]]] = ...,
+        point_scaling: str = ...,
+        title: str = ...,
+        save: T.Optional[str] = ...,
+    ) -> Figure: ...
+
+class Martini:
+    def __init__(
+        self,
+        source: T.Optional[SPHSource] = ...,
+        datacube: T.Optional[DataCube] = ...,
+        beam: T.Optional[_BaseBeam] = ...,
+        noise: T.Optional[_BaseNoise] = ...,
+        sph_kernel: T.Optional[_BaseSPHKernel] = ...,
+        spectral_model: T.Optional[_BaseSpectrum] = ...,
+        quiet: T.Optional[bool] = ...,
+    ) -> None: ...
+    @property
+    def datacube(self) -> DataCube: ...
     def insert_source_in_cube(
         self,
         skip_validation: bool = ...,
         progressbar: T.Optional[bool] = ...,
         ncpu: int = ...,
     ) -> None: ...
+    def convolve_beam(self) -> None: ...
+    def add_noise(self) -> None: ...
     def write_fits(
         self, filename: str, channels: str = ..., overwrite: bool = ...
     ) -> None: ...
@@ -51,14 +96,31 @@ class Martini:
         memmap: bool = ...,
         compact: bool = ...,
     ) -> None: ...
-    def reset(self) -> None: ...
-    def preview(
+
+class GlobalProfile(_BaseMartini):
+    def __init__(
+        self,
+        source: T.Optional[SPHSource] = ...,
+        spectral_model: T.Optional[_BaseSpectrum] = ...,
+        n_channels: int = ...,
+        channel_width: U.Quantity[U.km / U.s] = ...,
+        velocity_centre: U.Quantity[U.km / U.s] = ...,
+        channels: str = ...,
+        quiet: bool = ...,
+    ) -> None: ...
+    def insert_source_in_spectrum(self) -> None: ...
+    @property
+    def spectrum(self) -> U.Quantity[U.Jy]: ...
+    @property
+    def channel_edges(self) -> T.Union[U.Quantity[U.Hz], U.Quantity[U.km / U.s]]: ...
+    @property
+    def channel_mids(self) -> T.Union[U.Quantity[U.Hz], U.Quantity[U.km / U.s]]: ...
+    @property
+    def channel_width(self) -> U.Quantity[U.km / U.s]: ...
+    def plot_spectrum(
         self,
-        max_points: int = ...,
         fig: int = ...,
-        lim: T.Optional[T.Union[str, U.Quantity[U.kpc]]] = ...,
-        vlim: T.Optional[T.Union[str, U.Quantity[U.km / U.s]]] = ...,
-        point_scaling: str = ...,
         title: str = ...,
+        show_vsys: bool = ...,
         save: T.Optional[str] = ...,
     ) -> Figure: ...
diff --git a/martini/sources/sph_source.py b/martini/sources/sph_source.py
index 4b7cee8..a2d51d8 100644
--- a/martini/sources/sph_source.py
+++ b/martini/sources/sph_source.py
@@ -459,6 +459,7 @@ def preview(
         )
         size = 300 * size_scale if point_scaling == "auto" else 10
         fig = plt.figure(fig, figsize=(12, 4))
+        fig.clf()
         fig.suptitle(title)
 
         # ----- MOMENT 1 -----
diff --git a/martini/sph_kernels.py b/martini/sph_kernels.py
index 8f7953f..f988914 100644
--- a/martini/sph_kernels.py
+++ b/martini/sph_kernels.py
@@ -175,7 +175,8 @@ def _init_sm_ranges(self):
         """
         Determine maximum number of pixels reached by kernel.
         """
-        self.sm_ranges = np.ceil(self.sm_lengths * self.size_in_fwhm).astype(int)
+        # store as float: use case for np.inf in GlobalProfile:
+        self.sm_ranges = np.ceil(self.sm_lengths * self.size_in_fwhm)
 
         return
 
@@ -940,16 +941,22 @@ class DiracDeltaKernel(_BaseSPHKernel):
         0 &{\\rm for}\\;q > 0
         \\end{cases}
 
+    Parameters
+    ----------
+    size_in_fwhm : float
+        In principle the width of a DiracDelta kernel is 0, but this would
+        lead to no particles contributing to any pixels. Ideally this would
+        be set to approximately the size of the pixel, but setting it to
+        the smoothing length (1.0) is acceptable given the validation condition.
+        (Default: 1.0)
+
     """
 
     max_valid_size = 0.5
 
-    def __init__(self):
+    def __init__(self, size_in_fwhm=1.0):
         super().__init__()
-        # In principle the size for a DiracDelta kernel is 0, but this would
-        # lead to no particles being used. Ideally we would want ~the pixel
-        # size here, but the sph smoothing length is acceptable.
-        self.size_in_fwhm = 1
+        self.size_in_fwhm = size_in_fwhm
         self._rescale = 1  # need this to be present
         return
 
@@ -976,7 +983,7 @@ def kernel(self, q):
             Kernel value at positions q.
         """
 
-        return np.where(q, np.inf * np.ones(q.shape), np.zeros(q.shape))
+        return np.where(q, np.zeros(q.shape), np.inf * np.ones(q.shape))
 
     def _kernel_integral(self, dij, h, mask=np.s_[...]):
         """
@@ -998,7 +1005,6 @@ def _kernel_integral(self, dij, h, mask=np.s_[...]):
         out : array_like
             Kernel integral over the pixel area.
         """
-
         return np.where((np.abs(dij) < 0.5 * U.pix).all(axis=0), 1, 0) * U.pix**-2
 
     def _validate(self, sm_lengths, noraise=False, quiet=False):
diff --git a/martini/sph_kernels.pyi b/martini/sph_kernels.pyi
index 913ce30..0653efe 100644
--- a/martini/sph_kernels.pyi
+++ b/martini/sph_kernels.pyi
@@ -122,7 +122,7 @@ class DiracDeltaKernel(_BaseSPHKernel):
     max_valid_size: float
     size_in_fwhm: float
 
-    def __init__(self) -> None: ...
+    def __init__(self, size_in_fwhm: float = ...) -> None: ...
     def kernel(self, q: ndarray) -> ndarray: ...
     def _kernel_integral(
         self,


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index aa4e7e1..9b72312 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,7 +1,7 @@
 import pytest
 import numpy as np
 from astropy import units as U
-from martini.martini import Martini
+from martini.martini import Martini, GlobalProfile
 from martini.datacube import DataCube
 from martini.beams import GaussianBeam
 from martini.noise import GaussianNoise
@@ -337,3 +337,18 @@ def many_particle_source():
 @pytest.fixture(scope="function")
 def adaptive_kernel_test_source():
     yield adaptive_kernel_test_sourcegen
+
+
+@pytest.fixture(scope="function")
+def gp():
+    source = sps_sourcegen()
+    spectral_model = GaussianSpectrum()
+
+    m = GlobalProfile(
+        source=source,
+        spectral_model=spectral_model,
+        n_channels=16,
+        velocity_centre=source.distance * source.h * 100 * U.km / U.s / U.Mpc,
+    )
+    m.insert_source_in_spectrum()
+    yield m
diff --git a/tests/test_martini.py b/tests/test_martini.py
index b16dc20..856a78f 100644
--- a/tests/test_martini.py
+++ b/tests/test_martini.py
@@ -2,7 +2,7 @@
 import pytest
 import numpy as np
 import h5py
-from martini.martini import Martini
+from martini.martini import Martini, GlobalProfile, _BaseMartini
 from martini.datacube import DataCube, HIfreq
 from martini.beams import GaussianBeam
 from test_sph_kernels import simple_kernels
@@ -307,14 +307,32 @@ def test_add_noise(self, m_init):
             (-7 * U.km / U.s, False),
         ),
     )
+    @pytest.mark.parametrize("spatial", (True, False))
+    @pytest.mark.parametrize("spectral", (True, False))
     def test_prune_particles(
-        self, ra_off, ra_in, dec_off, dec_in, v_off, v_in, single_particle_source
+        self,
+        ra_off,
+        ra_in,
+        dec_off,
+        dec_in,
+        v_off,
+        v_in,
+        single_particle_source,
+        spatial,
+        spectral,
     ):
         """
         Check that a particle offset by a specific set of (RA, Dec, v) is inside/outside
         the cube as expected.
         """
-        expect_particle = all((ra_in, dec_in, v_in))
+        if spatial and spectral:
+            expect_particle = all((ra_in, dec_in, v_in))
+        elif spatial and not spectral:
+            expect_particle = all((ra_in, dec_in))
+        elif spectral and not spatial:
+            expect_particle = v_in
+        elif not spectral and not spatial:
+            expect_particle = True
         # set distance so that 1kpc = 1arcsec
         distance = (1 * U.kpc / 1 / U.arcsec).to(U.Mpc, U.dimensionless_angles())
         source = single_particle_source(
@@ -334,6 +352,7 @@ def test_prune_particles(
         beam = GaussianBeam(bmaj=1 * U.arcsec, bmin=1 * U.arcsec, truncate=4)
         sph_kernel = _CubicSplineKernel()
         spectral_model = GaussianSpectrum(sigma=1 * U.km / U.s)
+        # need to use _BaseMartini below to manipulate _prune_kwargs
         kwargs = dict(
             source=source,
             datacube=datacube,
@@ -341,6 +360,7 @@ def test_prune_particles(
             noise=None,
             sph_kernel=sph_kernel,
             spectral_model=spectral_model,
+            _prune_kwargs=dict(spatial=spatial, spectral=spectral),
         )
         # if more than 1px (datacube) + 5px (pad) + 2px (sm_range) then expect to prune
         # if more than 1px (datacube) + 4px (4*spectrum_half_width) then expect to prune
@@ -348,9 +368,9 @@ def test_prune_particles(
             with pytest.raises(
                 RuntimeError, match="No source particles in target region."
             ):
-                Martini(**kwargs)
+                _BaseMartini(**kwargs)
         else:
-            assert Martini(**kwargs).source.npart == 1
+            assert _BaseMartini(**kwargs).source.npart == 1
 
     def test_reset(self, m_nn):
         """
@@ -433,3 +453,174 @@ def test_parallel_consistent_with_serial(self, many_particle_source, dc_zeros):
         m.insert_source_in_cube(ncpu=2, progressbar=False)
 
         assert U.allclose(m.datacube._array, expected_result)
+
+
+class TestGlobalProfile:
+    @pytest.mark.parametrize("spectral_model", (DiracDeltaSpectrum, GaussianSpectrum))
+    @pytest.mark.parametrize("ra", (0 * U.deg, 180 * U.deg))
+    def test_mass_accuracy(self, spectral_model, single_particle_source, ra):
+        """
+        Check that the input mass in the particles gives the correct total mass in the
+        spectrum, by checking the conversion back to total mass. Covers testing
+        GlobalProfile.insert_source_in_spectrum.
+        """
+
+        # single_particle_source has a mass of 1E4Msun, temperature of 1E4K
+        # we test both ra=0deg and ra=180deg to make sure all particles included
+        source = single_particle_source(ra=ra)
+        m = GlobalProfile(
+            source=source,
+            spectral_model=spectral_model(),
+            n_channels=32,
+            channel_width=10 * U.km * U.s**-1,
+            velocity_centre=source.vsys,
+        )
+
+        m.insert_source_in_spectrum()
+
+        # flux
+        F = m.spectrum.sum()  # Jy
+
+        # distance
+        D = m.source.distance
+
+        # channel width
+        dv = m.channel_width
+
+        # HI mass
+        MHI = (
+            2.36e5
+            * U.Msun
+            * D.to_value(U.Mpc) ** 2
+            * F.to_value(U.Jy)
+            * dv.to_value(U.km / U.s)
+        ).to(U.Msun)
+
+        # demand accuracy within 1% after source insertion
+        assert U.isclose(MHI, m.source.mHI_g.sum(), rtol=1e-2)
+
+    @pytest.mark.parametrize(
+        ("ra_off", "ra_in"),
+        (
+            (0 * U.arcsec, True),
+            (3 * U.arcsec, True),
+            (5 * U.deg, False),  # global profile uses 1 deg pixel
+            (-3 * U.arcsec, True),
+            (-5 * U.deg, False),  # global profile uses 1 deg pixel
+        ),
+    )
+    @pytest.mark.parametrize(
+        ("dec_off", "dec_in"),
+        (
+            (0 * U.arcsec, True),
+            (3 * U.arcsec, True),
+            (5 * U.deg, False),  # global profile uses 1 deg pixel
+            (-3 * U.arcsec, True),
+            (-5 * U.deg, False),  # global profile uses 1 deg pixel
+        ),
+    )
+    @pytest.mark.parametrize(
+        ("v_off", "v_in"),
+        (
+            (0 * U.km / U.s, True),
+            (3 * U.km / U.s, True),
+            (7 * U.km / U.s, False),
+            (-3 * U.km / U.s, True),
+            (-7 * U.km / U.s, False),
+        ),
+    )
+    def test_prune_particles(
+        self, ra_off, ra_in, dec_off, dec_in, v_off, v_in, single_particle_source
+    ):
+        """
+        Check that a particle offset by a specific set of (RA, Dec, v) is inside/outside
+        the cube as expected. GlobalProfile should ignore RA, Dec when pruning.
+        """
+        # GlobalProfile should ignore RA, Dec when pruning:
+        expect_particle = v_in
+        # set distance so that 1kpc = 1arcsec
+        distance = (1 * U.kpc / 1 / U.arcsec).to(U.Mpc, U.dimensionless_angles())
+        source = single_particle_source(
+            distance=distance, ra=ra_off, dec=dec_off, vpeculiar=v_off
+        )
+        spectral_model = GaussianSpectrum(sigma=1 * U.km / U.s)
+        kwargs = dict(
+            source=source,
+            spectral_model=spectral_model,
+            n_channels=2,
+            channel_width=1 * U.km / U.s,
+            velocity_centre=source.distance * source.h * 100 * U.km / U.s / U.Mpc,
+        )
+        # if more than 1px (datacube) + 4px (4*spectrum_half_width) then expect to prune
+        if not expect_particle:
+            with pytest.raises(
+                RuntimeError, match="No source particles in target region."
+            ):
+                GlobalProfile(**kwargs)
+        else:
+            assert GlobalProfile(**kwargs).source.npart == 1
+
+    def test_reset(self, gp):
+        """
+        Check that resetting global profile instance zeros out datacube and spectrum.
+        """
+        cube_array = gp._datacube._array
+        assert gp._datacube._array.sum() > 0
+        spectrum = gp.spectrum
+        assert spectrum.sum() > 0
+        gp.reset()
+        assert gp._datacube._array.sum() == 0
+        assert not hasattr(gp, "_spectrum")
+        # check that can start over and get the same result w/o errors
+        gp.insert_source_in_spectrum()
+        assert U.allclose(cube_array, gp._datacube._array)
+        assert U.allclose(spectrum, gp.spectrum)
+        # check that can reset after doing nothing
+        gp.reset()
+        gp.reset()
+
+    def test_preview(self, gp):
+        """
+        Simply check that the preview visualisation runs without error.
+        """
+        # with default arguments
+        with pytest.warns(UserWarning, match="singular"):
+            # warning: single-particle source is used, so axis limits try to be equal
+            gp.preview()
+        # with non-default arguments
+        gp.preview(
+            max_points=1000,
+            fig=2,
+            lim="datacube",
+            vlim="datacube",
+            point_scaling="fixed",
+            title="test",
+        )
+
+    @pytest.mark.parametrize("channels", ("frequency", "velocity"))
+    def test_channel_modes(self, single_particle_source, channels):
+        """
+        Check that channels have expected units in both modes (frequency, velocity).
+        """
+        source = single_particle_source()
+        m = GlobalProfile(
+            source=source,
+            spectral_model=GaussianSpectrum(sigma="thermal"),
+            n_channels=32,
+            channel_width=10 * U.km * U.s**-1,
+            velocity_centre=source.vsys,
+            channels=channels,
+        )
+        expected_units = dict(frequency=U.Hz, velocity=U.km * U.s**-1)[channels]
+        # these will raise if there's a problem:
+        m.channel_edges.to(expected_units)
+        m.channel_mids.to(expected_units)
+
+    def test_view_spectrum(self, gp):
+        """
+        Simply check that plotting spectrum runs without error.
+        """
+        # with default arguments
+        gp.plot_spectrum()
+        # with non-default arguments
+        gp.plot_spectrum(fig=2, title="test", show_vsys=False)