diff --git a/docs/tutorials/spectrum.rst b/docs/tutorials/spectrum.rst index c054700..fba3c80 100644 --- a/docs/tutorials/spectrum.rst +++ b/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") diff --git a/pyproject.toml b/pyproject.toml index 083b7bf..26883d3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,7 @@ dependencies = [ "adjustText", "named-arrays~=2.8", "numba", + "scipy", "fiasco~=0.8.2", "sunpy~=8.0", ] diff --git a/utu/spectrum/__init__.py b/utu/spectrum/__init__.py index fe3d4b8..e5e9e56 100644 --- a/utu/spectrum/__init__.py +++ b/utu/spectrum/__init__.py @@ -7,6 +7,7 @@ ) from ._names import spectroscopic from ._plots import stem +from ._synthesis import synthesize __all__ = [ "contribution_function", @@ -14,4 +15,5 @@ "lines", "spectroscopic", "stem", + "synthesize", ] diff --git a/utu/spectrum/_lines.py b/utu/spectrum/_lines.py index abad7a3..7671a88 100644 --- a/utu/spectrum/_lines.py +++ b/utu/spectrum/_lines.py @@ -158,7 +158,7 @@ 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 @@ -166,6 +166,12 @@ def contribution_function( 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. diff --git a/utu/spectrum/_synthesis.py b/utu/spectrum/_synthesis.py new file mode 100644 index 0000000..5d733f5 --- /dev/null +++ b/utu/spectrum/_synthesis.py @@ -0,0 +1,411 @@ +"""The spectral radiance of an optically thin plasma.""" + +from typing import cast + +import astropy.constants +import astropy.units as u +import named_arrays as na +import numpy as np +import scipy.special + +__all__ = [ + "synthesize", +] + +_size_block = 2**25 +""" +The most elements :func:`synthesize` puts in any one array along the way, a +quarter of a gigabyte in double precision. + +The largest of them is the product of the contribution function, the +emission measure, and the profile of a block of lines in one bin of +velocity, at every temperature and in every pixel. The lines are synthesized +as many at a time as fit in this, and at least one. +""" + + +def synthesize( + contribution: na.FunctionArray, + temperature: na.AbstractScalar, + mass: u.Quantity | na.AbstractScalar, + emission_measure: u.Quantity | na.AbstractScalar, + velocity: na.AbstractScalar, + shift: u.Quantity | na.AbstractScalar = 0 * u.km / u.s, + width_nonthermal: u.Quantity | na.AbstractScalar = 0 * u.km / u.s, + axis_temperature: str = "temperature", + axis_velocity: str = "velocity", + axis: str = "line", +) -> na.FunctionArray[na.DopplerVectorArray, na.AbstractScalar]: + r""" + Compute the spectral radiance of an optically thin plasma, line by line. + + Each line is spread over a range of velocities around its rest + wavelength, as a Gaussian whose center is the line-of-sight velocity of + the plasma and whose width is that of the thermal motion of its ions + added in quadrature to a nonthermal width. The plasma at each temperature + has a profile of its own, and the spectrum is their sum, + + .. math:: + + L(\lambda) = \frac{1}{4 \pi} \sum_T G(T) \, \mathrm{EM}(T) \, + \phi(\lambda; T), + + where :math:`G` is the contribution function, :math:`\mathrm{EM}` the + emission measure, and :math:`\phi` the profile, normalized to unit area. + Its standard deviation, in velocity, is + + .. math:: + + \sigma^2 = \frac{k_B T}{m} + \sigma_\mathrm{nt}^2, + + with :math:`m` the mass of the ion and :math:`\sigma_\mathrm{nt}` the + nonthermal width. + + The radiance of each bin of velocity is the integral of the profile + across it divided by its width in wavelength, not the profile sampled at + its center. So the sum of the radiance over the bins is the radiance of + the line, however coarse the bins, as long as they span the line. The + integral is taken from whichever tail of the profile the bin is in, so + it keeps its precision far into the wings. + + Parameters + ---------- + contribution + The contribution functions of the lines, as + :func:`~utu.spectrum.contribution_function` returns them: the rest + wavelength of each line along `axis` as the inputs, and its + contribution function along `axis_temperature` and `axis` as the + outputs. :mod:`fiasco` leaves the ratio of hydrogen to electrons out + of these, so the emission measure they go with is that of + :math:`n_e n_H`. + temperature + The temperatures of the plasma, along `axis_temperature`: the ones + `contribution` was computed at. They must be positive. + mass + The mass of the ion each line belongs to, such as + :attr:`fiasco.Ion.mass`. Give it along `axis` to synthesize the lines + of several ions at once. + emission_measure + How much plasma there is in each bin of temperature, + :math:`\int n_e n_H \, dh` through the column, along + `axis_temperature`. Every other axis is carried through to the + result, such as the two axes of a map. + + This is the emission measure of each bin, not the differential + emission measure: a DEM per unit :math:`\log_{10} T`, such as + :func:`utu.dem.plowman` returns, has the same units, and has to be + multiplied by the width of each bin in :math:`\log_{10} T` first. + velocity + The vertices of the bins of velocity, along `axis_velocity`. Positive + velocities are away from the observer, to longer wavelengths. + shift + The line-of-sight velocity of the plasma. It may vary along any axis + of `emission_measure`, including `axis_temperature`, and along + `axis`. + width_nonthermal + The standard deviation of the nonthermal motions of the plasma along + the line of sight. Often quoted instead as the most probable speed + :math:`\xi`, which is :math:`\sqrt{2}` times larger. It may vary + along the same axes as `shift`. + axis_temperature + The name of the axis of `temperature`, summed over. + axis_velocity + The name of the axis of `velocity`. + axis + The name of the axis along the lines of `contribution`. Any other + argument which varies along it must have one element along it, or as + many as `contribution`. + + Returns + ------- + : + The spectral radiance of each line in each bin of velocity, along + `axis`, `axis_velocity`, and every axis of the other arguments but + `axis_temperature`. The inputs are the vertices of the bins, as + wavelengths with the rest wavelength of each line beside them. + + Raises + ------ + ValueError + If `velocity` has no axis `axis_velocity` or varies along + `axis_temperature`, if `temperature` or the + outputs of `contribution` have no axis `axis_temperature`, if a + temperature is not positive, or if an argument varies along `axis` + differently than `contribution` does. + + Notes + ----- + Nothing here is particular to an instrument: the profile is that of the + plasma alone, with no instrumental width, and every line in + `contribution` is synthesized, so choosing which lines are worth + including is left to the caller. Natural and pressure broadening are + left out, which in the transition region and the corona is a small + fraction of a milliangstrom, and so is opacity. + + The result is filled in one bin of velocity and one block of lines at a + time, with as many lines in a block as keep the largest array along the + way under :data:`_size_block` elements. So a map of many lines needs + little more memory than the result itself. Where the profile is the same + for every line, because none of `mass`, `shift`, and `width_nonthermal` + varies along `axis`, it is computed once for all of them. + + Examples + -------- + The :math:`\mathrm{Ne\,VII}\;465\,\AA` line of plasma near half a million + kelvin, at rest and moving toward the observer at + :math:`20\,\mathrm{km\,s^{-1}}`. + + .. jupyter-execute:: + + import astropy.units as u + import fiasco + import matplotlib.pyplot as plt + import named_arrays as na + import numpy as np + import utu + + fiasco.log.setLevel("ERROR") + + axis = "temperature" + temperature = na.geomspace(1e5, 1e6, axis=axis, num=21) * u.K + + # a Gaussian in the logarithm of the temperature, holding + # 1e27 cm^-5 of plasma in all + weight = np.exp(-np.square((np.log10(temperature / u.K) - 5.7) / 0.1) / 2) + emission_measure = 1e27 / u.cm**5 * weight / weight.sum(axis) + + ion = fiasco.Ion("Ne 7", temperature.ndarray) + + contribution = utu.spectrum.contribution_function( + ion=ion, + density=1e15 * u.K / u.cm**3 / temperature, + axis_temperature=axis, + ) + + # the strongest line of the ion + brightest = np.argmax(contribution.outputs.max(axis), axis="line") + contribution = contribution[brightest] + + shift = na.ScalarArray([0, -20] * u.km / u.s, axes="shift") + + result = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=ion.mass, + emission_measure=emission_measure, + velocity=na.linspace(-100, 100, axis="velocity", num=81) * u.km / u.s, + shift=shift, + axis_temperature=axis, + ) + + fig, ax = plt.subplots(figsize=(6, 4), constrained_layout=True) + na.plt.stairs( + result.inputs.velocity, + result.outputs, + ax=ax, + axis="velocity", + label=shift, + ) + ax.set_xlabel(f"velocity ({u.km / u.s:latex_inline})") + ax.set_ylabel(f"spectral radiance\n({result.outputs.unit:latex_inline})") + ax.legend(title="shift"); + """ + temperature = _scalar(temperature) + mass = _scalar(mass) + emission_measure = _scalar(emission_measure) + velocity = _scalar(velocity).explicit + shift = _scalar(shift) + width_nonthermal = _scalar(width_nonthermal) + + g = contribution.outputs + wavelength_rest = contribution.inputs + + if axis_velocity not in na.shape(velocity): + raise ValueError( + f"`velocity` must have the axis {axis_velocity!r}, " + f"got {na.shape(velocity)}." + ) + if axis_temperature in na.shape(velocity): + raise ValueError( + f"`velocity` must not vary along {axis_temperature!r}, which is " + f"summed over, got {na.shape(velocity)}." + ) + for name, a in dict(temperature=temperature, contribution=g).items(): + if axis_temperature not in na.shape(a): + raise ValueError( + f"`{name}` must have the axis {axis_temperature!r}, " + f"got {na.shape(a)}." + ) + if not np.all(temperature > 0 * u.K): + raise ValueError("every temperature must be positive.") + + # The number of lines, or None for a single line without an axis of + # lines. Every other argument has to agree with it, or a line would be + # profiled with the mass or the motions of another. + num_line = na.shape(g).get(axis) + profile = dict( + mass=mass, + velocity=velocity, + shift=shift, + width_nonthermal=width_nonthermal, + ) + for name, a in (profile | dict(emission_measure=emission_measure)).items(): + num = na.shape(a).get(axis) + if num is None: + continue + if num_line is None: + raise ValueError( + f"`{name}` varies along {axis!r}, " + f"but `contribution` is a single line." + ) + if num not in (1, num_line): + raise ValueError( + f"`{name}` has {num} elements along {axis!r}, " + f"but `contribution` has {num_line} lines." + ) + + shape = na.shape_broadcasted( + g, + wavelength_rest, + temperature, + emission_measure, + *profile.values(), + ) + num_temperature = shape.pop(axis_temperature) + num_bin = shape[axis_velocity] - 1 + shape[axis_velocity] = num_bin + + unit = u.erg / u.cm**2 / u.s / u.sr / u.AA + result = na.ScalarArray( + ndarray=np.zeros(tuple(shape.values())) << unit, + axes=tuple(shape), + ) + + # Enough lines at a time that the product of the contribution function, + # the emission measure, and the profile in one bin, at every temperature + # and in every pixel, stays under the size of a block. + size_line = num_temperature + for a, n in shape.items(): + if a not in (axis, axis_velocity): + size_line *= n + num_block = max(1, _size_block // size_line) + if num_line is None: + blocks = [None] + else: + blocks = [slice(i, i + num_block) for i in range(0, num_line, num_block)] + + # The profile is the same for every line unless something it depends on + # varies along the axis of lines, and is then computed once per bin. + shared = all(axis not in na.shape(a) for a in profile.values()) + + for k in range(num_bin): + index_bin = {axis_velocity: slice(k, k + 2)} + + if shared: + fraction_shared = _fraction( + temperature=temperature, + **{name: a[index_bin] for name, a in profile.items()}, + axis_velocity=axis_velocity, + ) + + for block in blocks: + index = {} if block is None else {axis: block} + + if shared: + fraction = fraction_shared + else: + fraction = _fraction( + temperature=temperature, + **{ + name: _select(a[index_bin], index) + for name, a in profile.items() + }, + axis_velocity=axis_velocity, + ) + + radiance = g[index] * _select(emission_measure, index) * fraction + radiance = radiance.sum(axis_temperature) + + # The width of the bin in wavelength, which turns the radiance in + # it into a spectral radiance. The Doppler shift is linear in + # velocity, so this is exact, and does not take the difference of + # two wavelengths which agree in most of their digits. + width = np.diff(_select(velocity[index_bin], index), axis=axis_velocity) + width = wavelength_rest[index] * width / astropy.constants.c + + radiance = radiance / width / (4 * np.pi * u.sr) + + result[index | {axis_velocity: slice(k, k + 1)}] = radiance + + return na.FunctionArray( + inputs=na.DopplerVectorArray.from_velocity( + velocity=velocity, + wavelength_rest=wavelength_rest, + ), + outputs=result, + ) + + +def _fraction( + temperature: na.AbstractScalar, + mass: na.AbstractScalar, + velocity: na.AbstractScalar, + shift: na.AbstractScalar, + width_nonthermal: na.AbstractScalar, + axis_velocity: str, +) -> na.AbstractScalar: + """ + The fraction of the profile of a line which falls in one bin of velocity. + + `velocity` is the two vertices of the bin, along `axis_velocity`, and the + result keeps that axis, with one element. + + The fraction is the difference of the cumulative distribution of the + profile at the two vertices. Far into either wing that distribution is + within rounding of zero or of one, and the difference of two values near + one loses every digit, so it is taken instead from the complementary + error function of the tail the bin is in, which keeps its relative + precision there. + """ + width = np.sqrt( + astropy.constants.k_B * temperature / mass + np.square(width_nonthermal) + ) + + z = (velocity - shift) / (np.sqrt(2) * width) + z = z.to(u.dimensionless_unscaled).value + + z0 = z[{axis_velocity: slice(0, 1)}] + z1 = z[{axis_velocity: slice(1, 2)}] + + # +1 if the middle of the bin is above the center of the profile, -1 if + # below, so that the error function is evaluated in its own tail + side = np.where(z0 + z1 >= 0, 1, -1) + + erfc0 = scipy.special.erfc(side * z0) + erfc1 = scipy.special.erfc(side * z1) + + return side * (erfc0 - erfc1) / 2 + + +def _select( + a: na.AbstractScalar, + index: dict[str, slice], +) -> na.AbstractScalar: + """ + Select a block of lines from an argument of :func:`synthesize`, unless + it broadcasts along the axis of lines. + + An argument with one element along that axis is the same for every + line, and indexing it with a block past its first element would fail. + """ + for axis in index: + if na.shape(a).get(axis, 1) == 1: + return a + return a[index] + + +def _scalar( + a: u.Quantity | na.AbstractScalar, +) -> na.AbstractScalar: + """A quantity as a named scalar, so that it can be indexed by axis.""" + return cast(na.AbstractScalar, na.as_named_array(a)) diff --git a/utu/spectrum/_tests/test_synthesis.py b/utu/spectrum/_tests/test_synthesis.py new file mode 100644 index 0000000..5746a59 --- /dev/null +++ b/utu/spectrum/_tests/test_synthesis.py @@ -0,0 +1,399 @@ +import astropy.constants +import astropy.units as u +import fiasco +import named_arrays as na +import numpy as np +import pytest +import scipy.special + +import utu +import utu.spectrum._synthesis + +from .test_lines import needs_database + +# Two lines of two ions, with contribution functions made up rather than read +# from the database, so that these tests run wherever the package does. The +# shape of a contribution function has nothing to do with how it is spread +# into a profile. +temperature = na.geomspace(1e5, 10**6.5, axis="temperature", num=6) * u.K + +contribution = na.FunctionArray( + inputs=na.ScalarArray([465.221, 499.406] * u.AA, axes="line"), + outputs=na.ScalarArray( + ndarray=np.array( + [ + [1, 3, 2, 1, 0.5, 0.1], + [0.1, 0.2, 1, 4, 3, 1], + ] + ) + * 1e-24 + * u.erg + * u.cm**3 + / u.s, + axes=("line", "temperature"), + ), +) + +mass = na.ScalarArray([20.18, 28.09] * u.u, axes="line") + +emission_measure = na.ScalarArray( + ndarray=np.array([1, 2, 4, 2, 1, 0.5]) * 1e26 / u.cm**5, + axes="temperature", +) + +# wide enough to hold eleven standard deviations of the broadest line here, +# Ne VII at the hottest temperature, so that nothing is lost off the ends +velocity = na.linspace(-400, 400, axis="velocity", num=801) * u.km / u.s + + +def _synthesize(**kwargs: object) -> na.FunctionArray: + """:func:`utu.spectrum.synthesize` of the lines above, unless overridden.""" + arguments = dict( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=emission_measure, + velocity=velocity, + ) + return utu.spectrum.synthesize(**(arguments | kwargs)) + + +def _radiance(result: na.FunctionArray, axis: str = "velocity") -> na.AbstractScalar: + """The radiance in each bin of velocity, rather than per unit wavelength.""" + return result.outputs * np.diff(result.inputs.wavelength, axis=axis) + + +def _total(result: na.FunctionArray) -> na.AbstractScalar: + """The radiance of each line, summed over the bins of velocity.""" + return _radiance(result).sum("velocity") + + +def _moments(result: na.FunctionArray) -> tuple[na.AbstractScalar, na.AbstractScalar]: + """The mean and the variance of the velocity of each line.""" + v = result.inputs.velocity + v = (v[{"velocity": slice(1, None)}] + v[{"velocity": slice(None, -1)}]) / 2 + weight = result.outputs / result.outputs.sum("velocity") + mean = (weight * v).sum("velocity") + variance = (weight * np.square(v - mean)).sum("velocity") + return mean, variance + + +def test_synthesize(): + result = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=emission_measure, + velocity=velocity, + ) + + assert isinstance(result, na.FunctionArray) + assert isinstance(result.inputs, na.DopplerVectorArray) + assert result.outputs.unit.is_equivalent(u.erg / u.cm**2 / u.s / u.sr / u.AA) + assert np.all(result.outputs >= 0) + + # every line has its own bins, at the velocities asked for + assert na.shape(result.outputs) == dict(line=2, velocity=800) + assert np.allclose(result.inputs.velocity, velocity) + assert np.all(result.inputs.wavelength_rest == contribution.inputs) + + +def test_synthesize_conserves(): + """The bins of a line hold all of its radiance, however coarse they are.""" + expected = (contribution.outputs * emission_measure).sum("temperature") + expected = expected / (4 * np.pi * u.sr) + + for num in [801, 9]: + result = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=emission_measure, + velocity=na.linspace(-400, 400, axis="velocity", num=num) * u.km / u.s, + ) + assert np.allclose(_total(result), expected, rtol=1e-10) + + +def test_synthesize_symmetric(): + """A line at rest is the same on either side of its rest wavelength.""" + result = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=emission_measure, + velocity=velocity, + ) + + reverse = result.outputs[{"velocity": slice(None, None, -1)}] + assert np.allclose(reverse, result.outputs, rtol=1e-10) + + +@pytest.mark.parametrize( + argnames="shift", + argvalues=[-30, 0, 45] * u.km / u.s, +) +def test_synthesize_shift(shift: u.Quantity): + result = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=emission_measure, + velocity=velocity, + shift=shift, + ) + + mean, _ = _moments(result) + assert np.allclose(mean, shift, atol=1e-6 * u.km / u.s) + + +@pytest.mark.parametrize( + argnames="width_nonthermal", + argvalues=[0, 20] * u.km / u.s, +) +def test_synthesize_width(width_nonthermal: u.Quantity): + """Plasma at one temperature has the width of its ions' thermal motion.""" + index = {"temperature": slice(3, 4)} + + result = utu.spectrum.synthesize( + contribution=contribution[index], + temperature=temperature[index], + mass=mass, + emission_measure=emission_measure[index], + velocity=velocity, + width_nonthermal=width_nonthermal, + ) + + _, variance = _moments(result) + + expected = astropy.constants.k_B * temperature[index] / mass + expected = expected + np.square(width_nonthermal) + + # Sheppard's correction: binning a Gaussian adds the variance of a uniform + # distribution across one bin + step = velocity[{"velocity": 1}] - velocity[{"velocity": 0}] + expected = expected + np.square(step) / 12 + + assert np.allclose(variance, expected, rtol=1e-6) + + +def test_synthesize_single(): + """A line on its own, without an axis of lines, is the same line.""" + result = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=emission_measure, + velocity=velocity, + ) + + for i in range(2): + index = {"line": i} + single = utu.spectrum.synthesize( + contribution=contribution[index], + temperature=temperature, + mass=mass[index], + emission_measure=emission_measure, + velocity=velocity, + ) + assert "line" not in na.shape(single.outputs) + assert np.allclose(single.outputs, result.outputs[index], rtol=1e-12) + + +def test_synthesize_map(): + """ + The axes of the emission measure are carried through, and the shift and + the nonthermal width can vary along them and with temperature, each pixel + on its own. + """ + em = emission_measure * na.ScalarArray(np.array([1, 2, 3]), axes="x") + em = em * na.ScalarArray(np.array([1, 0.5]), axes="y") + + shift = na.ScalarArray([-20, 0, 30] * u.km / u.s, axes="x") + shift = shift + na.linspace(0, 10, axis="temperature", num=6) * u.km / u.s + + width = na.ScalarArray([5, 25] * u.km / u.s, axes="y") + width = width + na.linspace(0, 10, axis="temperature", num=6) * u.km / u.s + + result = _synthesize( + emission_measure=em, + shift=shift, + width_nonthermal=width, + ) + + assert na.shape(result.outputs) == dict(line=2, velocity=800, x=3, y=2) + + for i in range(3): + for j in range(2): + index = {"x": i, "y": j} + pixel = _synthesize( + emission_measure=em[index], + shift=shift[index], + width_nonthermal=width[index], + ) + assert np.allclose(pixel.outputs, result.outputs[index], rtol=1e-12) + + +def test_synthesize_axes(): + """Axes named other than the defaults give the same result.""" + rename = dict(temperature="t", velocity="v", line="l") + + def _rename(a: na.AbstractScalar) -> na.AbstractScalar: + """The same array, with every axis renamed.""" + a = a.explicit + return na.ScalarArray(a.ndarray, axes=tuple(rename[x] for x in a.axes)) + + result = utu.spectrum.synthesize( + contribution=na.FunctionArray( + inputs=_rename(contribution.inputs), + outputs=_rename(contribution.outputs), + ), + temperature=_rename(temperature), + mass=_rename(mass), + emission_measure=_rename(emission_measure), + velocity=_rename(velocity), + axis_temperature="t", + axis_velocity="v", + axis="l", + ) + + expected = _synthesize() + + assert na.shape(result.outputs) == dict(l=2, v=800) + assert np.allclose(_rename(expected.outputs), result.outputs, rtol=1e-12) + + +@pytest.mark.parametrize( + argnames="kwargs", + argvalues=[ + # the vertices of the bins along some other axis than the default + dict(velocity=na.linspace(-400, 400, axis="v", num=801) * u.km / u.s), + # bins which change with the temperature summed over + dict(velocity=velocity * na.linspace(1, 2, axis="temperature", num=6)), + # temperatures along some other axis than the default + dict(temperature=na.geomspace(1e5, 1e6, axis="t", num=6) * u.K), + # plasma with no temperature + dict(temperature=na.linspace(0, 1e6, axis="temperature", num=6) * u.K), + # more masses than lines + dict(mass=na.ScalarArray([20, 28, 32] * u.u, axes="line")), + # a mass for each of two lines, given to one line + dict(contribution=contribution[{"line": 0}]), + ], + ids=["velocity", "velocity-temperature", "temperature", "zero", "mass", "single"], +) +def test_synthesize_error(kwargs: dict): + with pytest.raises(ValueError): + _synthesize(**kwargs) + + +def test_synthesize_empty(): + """No lines at all, such as a passband with none in it, is no lines.""" + index = {"line": slice(0, 0)} + + result = _synthesize( + contribution=contribution[index], + mass=mass[index], + ) + + assert na.shape(result.outputs) == dict(line=0, velocity=800) + + +def test_synthesize_broadcast(): + """A mass with one element along the axis of lines is every line's.""" + one = na.ScalarArray([20.18] * u.u, axes="line") + + result = _synthesize(mass=one) + expected = _synthesize(mass=20.18 * u.u) + + assert np.allclose(result.outputs, expected.outputs, rtol=1e-12) + + +def test_synthesize_reversed(): + """Bins from fast to slow are the same bins as from slow to fast.""" + reverse = {"velocity": slice(None, None, -1)} + + result = _synthesize(velocity=velocity.explicit[reverse]) + expected = _synthesize() + + assert np.all(result.outputs > 0) + assert np.allclose(result.outputs[reverse], expected.outputs, rtol=1e-10) + + +def test_synthesize_wings(): + """ + The far wings of a line hold what the tail of a Gaussian holds, rather + than nothing at all. + """ + index = {"temperature": slice(3, 4)} + sigma = np.sqrt(astropy.constants.k_B * temperature[{"temperature": 3}] / mass) + sigma = sigma.to(u.km / u.s) + + # vertices every tenth of a standard deviation, out to thirty of them + num = 601 + z = na.linspace(-30, 30, axis="velocity", num=num) + + result = _synthesize( + contribution=contribution[index], + temperature=temperature[index], + emission_measure=emission_measure[index], + velocity=z * sigma, + ) + + radiance = _radiance(result) + assert np.all(radiance[{"velocity": slice(50, -50)}] > 0) + + # the fraction of the line beyond eight standard deviations, about 6e-16 + wing = {"velocity": slice((num - 1) // 2 + 80, None)} + fraction = radiance[wing].sum("velocity") / radiance.sum("velocity") + expected = scipy.special.erfc(8 / np.sqrt(2)) / 2 + assert np.allclose(fraction, expected, rtol=1e-6) + + +def test_synthesize_blocks(monkeypatch: pytest.MonkeyPatch): + """A map synthesized a line at a time is the same map.""" + em = emission_measure * na.ScalarArray(np.array([1, 2, 3]), axes="x") + expected = _synthesize(emission_measure=em) + + # room for one line at every temperature and in every pixel, and no more + monkeypatch.setattr(utu.spectrum._synthesis, "_size_block", 6 * 3) + result = _synthesize(emission_measure=em) + + assert np.allclose(result.outputs, expected.outputs, rtol=1e-12) + + +@needs_database +def test_synthesize_lines(): + """The radiance of the lines is that :func:`utu.spectrum.lines` gives.""" + t = na.geomspace(1e5, 1e6, axis="temperature", num=6) * u.K + density = 1e15 * u.K / u.cm**3 / t + em = 1e27 / u.cm**5 * na.ScalarArray(np.ones(6), axes="temperature") + + ratio = fiasco.proton_electron_ratio(t.ndarray) + + ion = fiasco.Ion("O 5", t.ndarray) + + g = utu.spectrum.contribution_function( + ion=ion, + density=density, + axis_temperature="temperature", + proton_electron_ratio=ratio, + ) + + result = utu.spectrum.synthesize( + contribution=g, + temperature=t, + mass=ion.mass, + emission_measure=em, + velocity=velocity, + ) + + expected = utu.spectrum.lines( + temperature=t, + density=density, + emission_measure=em, + ions=["O 5"], + proton_electron_ratio=ratio, + ) + + total = _total(result).sum("line") + expected = expected.outputs.sum("line") / (4 * np.pi * u.sr) + assert np.allclose(total, expected, rtol=1e-10)