Add utu.spectrum.synthesize, the spectral radiance of each line

docs/tutorials/spectrum.rst

@@ -59,7 +59,10 @@ density per temperature, and a density on an axis of its own as every density
 at every temperature.
 
 The emission measure is how much material there is at each temperature,
-:math:`\int n_e^2 \, dh` through the column. A real one is measured; this one
+:math:`\int n_e n_H \, dh` through the column. It is the product of the
+densities of electrons and of hydrogen rather than the square of the first,
+since :mod:`fiasco` leaves the ratio of the two out of its contribution
+functions. A real one is measured; this one
 is a Gaussian in :math:`\log T` centred where the transition region emits,
 holding :math:`10^{27}\,\mathrm{cm}^{-5}` of material in all.
 
@@ -185,3 +188,72 @@ radiance computed above has to be put back into one before it can be drawn.
     )
     ax.set_xlabel(f"wavelength ({u.AA:latex_inline})")
     ax.set_ylabel(f"radiance ({radiance.unit:latex_inline})");
+
+
+Line profiles
+-------------
+
+A line is not infinitely narrow. The ions emitting it move, so each photon is
+Doppler shifted by the velocity of its ion along the line of sight, and the
+line is spread over a range of wavelengths about as wide as the motions of
+the plasma: the thermal motion of its ions, any unresolved nonthermal motion,
+and a shift of the whole profile if the plasma is flowing.
+:func:`~utu.spectrum.synthesize` spreads each line into a Gaussian with
+these properties, in bins of velocity around its rest wavelength. The plasma
+at each temperature gets a profile of its own, as wide as the thermal motion
+at that temperature, so the hot and the cool parts of a column can move
+differently.
+
+It takes the contribution functions of an ion as
+:func:`~utu.spectrum.contribution_function` returns them, along with the
+mass of the ion, which sets its thermal width. Here is the
+:math:`\mathrm{O\,V}` line from above, from plasma with no nonthermal
+motion and from plasma with a nonthermal width of
+:math:`20\,\mathrm{km\,s^{-1}}`.
+
+.. jupyter-execute::
+
+    ion = fiasco.Ion("O 5", temperature.ndarray)
+
+    contribution = utu.spectrum.contribution_function(
+        ion=ion,
+        density=density,
+        axis_temperature=axis,
+    )
+
+    index = np.argmax(contribution.outputs.max(axis), axis="line")
+
+    width_nonthermal = na.ScalarArray([0, 20] * u.km / u.s, axes="width")
+
+    result = utu.spectrum.synthesize(
+        contribution=contribution[index],
+        temperature=temperature,
+        mass=ion.mass,
+        emission_measure=emission_measure,
+        velocity=na.linspace(-150, 150, axis="velocity", num=121) * u.km / u.s,
+        width_nonthermal=width_nonthermal,
+        axis_temperature=axis,
+    )
+
+    fig, ax = plt.subplots(figsize=(7, 3.5), constrained_layout=True)
+    na.plt.stairs(
+        result.inputs.wavelength,
+        result.outputs,
+        ax=ax,
+        axis="velocity",
+        label=width_nonthermal,
+    )
+    ax.set_xlabel(f"wavelength ({u.AA:latex_inline})")
+    ax.set_ylabel(f"spectral radiance\n({result.outputs.unit:latex_inline})")
+    ax.legend(title="nonthermal width");
+
+The radiance of each bin is the integral of the profile across it, so
+broadening a line moves its light around without adding to it, and loses
+none of it as long as the bins span the line. Summed over the bins, both
+profiles give the radiance of the line in the table above.
+
+.. jupyter-execute::
+
+    width_bin = np.diff(result.inputs.wavelength, axis="velocity")
+
+    (result.outputs * width_bin).sum("velocity")

pyproject.toml

@@ -22,6 +22,7 @@ dependencies = [
     "adjustText",
     "named-arrays~=2.8",
     "numba",
+    "scipy",
     "fiasco~=0.8.2",
     "sunpy~=8.0",
 ]

utu/spectrum/__init__.py

@@ -7,11 +7,13 @@
 )
 from ._names import spectroscopic
 from ._plots import stem
+from ._synthesis import synthesize
 
 __all__ = [
     "contribution_function",
     "ions",
     "lines",
     "spectroscopic",
     "stem",
+    "synthesize",
 ]

utu/spectrum/_lines.py

@@ -158,14 +158,20 @@ def contribution_function(
     axis: str = "line",
     proton_electron_ratio: None | u.Quantity | na.AbstractScalar = None,
 ) -> na.FunctionArray:
-    """
+    r"""
     Compute the contribution function of every line of an ion.
 
     Returned as a function of wavelength, so that a line and its strength
     cannot come apart. They are separate arrays underneath, of different
     lengths whenever an ion has a two-photon transition, and pairing them by
     hand is a way to label every line with its neighbor's wavelength.
 
+    The contribution function is that of :mod:`fiasco`, which leaves out the
+    ratio of the density of hydrogen to that of electrons. So the emission
+    measure it is multiplied by to find an intensity is
+    :math:`\int n_e n_H \, dh`, not :math:`\int n_e^2 \, dh`, which would
+    overstate the intensity by :math:`n_e / n_H`, about 1.2 in the corona.
+
     Whether the temperatures and the densities are taken in pairs or as a
     grid is decided by their axes. A density which shares the axis of the
     temperature describes one density per temperature, an isobaric
@@ -272,7 +278,7 @@ def lines(
     axis: str = "line",
     **kwargs: object,
 ) -> na.FunctionArray:
-    """
+    r"""
     Compute the emission lines of an optically thin plasma, brightest first.
 
     Every line of every ion abundant enough to contribute, with the intensity
@@ -298,7 +304,12 @@ def lines(
         atmosphere for instance; a density on its own axis is every density
         at every temperature.
     emission_measure
-        How much plasma there is at each temperature.
+        How much plasma there is at each temperature,
+        :math:`\int n_e n_H \, dh` through the column in each bin of
+        temperature. It is the product of the densities of electrons and of
+        hydrogen rather than the square of the first, since :mod:`fiasco`
+        leaves the ratio of the two out of the contribution functions, see
+        :func:`contribution_function`.
     wavelength_min
         The shortest wavelength worth computing.
         If :obj:`None` (the default), there is no lower bound.