From 635aabf53a2f632f80073cf39c13e2af19f7f541 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Wed, 14 Oct 2026 20:16:21 -0600 Subject: [PATCH 1/2] Add `utu.spectrum.synthesize`, the spectral radiance of each line Spread each line of `contribution_function()` into a Gaussian in bins of velocity around its rest wavelength. Its width is the thermal motion of the ion added in quadrature to a nonthermal width, and its center is the line-of-sight velocity of the plasma. The plasma at each temperature has a profile of its own, and the profiles are summed over temperature with the emission measure. Every other axis of the emission measure, such as the two axes of a map, is carried through, and the shift and the nonthermal width can vary along any of them. The radiance of each bin is the integral of the profile across it, so the radiance of a line is conserved on coarse grids. The result's inputs are a `DopplerVectorArray` of the bin vertices. The tutorial gains a section on line profiles, and now gives the emission measure as that of n_e n_H, which is what fiasco's contribution function goes with. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_017xqpht97bF4FWsACAA4wQa --- docs/tutorials/spectrum.rst | 78 ++++++- pyproject.toml | 1 + utu/spectrum/__init__.py | 2 + utu/spectrum/_synthesis.py | 279 ++++++++++++++++++++++++++ utu/spectrum/_tests/test_synthesis.py | 257 ++++++++++++++++++++++++ 5 files changed, 616 insertions(+), 1 deletion(-) create mode 100644 utu/spectrum/_synthesis.py create mode 100644 utu/spectrum/_tests/test_synthesis.py diff --git a/docs/tutorials/spectrum.rst b/docs/tutorials/spectrum.rst index c054700..78dff47 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,76 @@ 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:: + + wavelength = result.inputs.wavelength + width_bin = ( + wavelength[{"velocity": slice(1, None)}] + - wavelength[{"velocity": slice(None, -1)}] + ) + + (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/_synthesis.py b/utu/spectrum/_synthesis.py new file mode 100644 index 0000000..26e01e3 --- /dev/null +++ b/utu/spectrum/_synthesis.py @@ -0,0 +1,279 @@ +"""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", +] + + +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. + + 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. + 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 at each temperature, + :math:`\int n_e n_H \, dh` through the column in each bin of + temperature, along `axis_temperature`. Every other axis is carried + through to the result, such as the two axes of a map. + 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`. + + Returns + ------- + : + The spectral radiance of each line in each bin of velocity, along + `axis`, `axis_velocity`, and every axis of `emission_measure` but + `axis_temperature`. The inputs are the vertices of the bins, as + wavelengths with the rest wavelength of each line beside them. + + 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 lines are synthesized one at a time, so that the largest array held + at once is the profile of one line at every temperature and in every + pixel, rather than that of every line. + + 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"); + """ + wavelength_rest = contribution.inputs + + arrays = dict( + temperature=_scalar(temperature), + mass=_scalar(mass), + emission_measure=_scalar(emission_measure), + velocity=_scalar(velocity), + shift=_scalar(shift), + width_nonthermal=_scalar(width_nonthermal), + ) + + shape = na.shape(contribution.outputs) + + # A single line, such as one picked out of the lines of an ion with an + # integer index, has no axis of lines to loop over. + if axis not in shape: + result = _line( + contribution=contribution, + **arrays, + axis_temperature=axis_temperature, + axis_velocity=axis_velocity, + ) + else: + result = na.stack( + [ + _line( + contribution=contribution[{axis: i}], + **{name: a[{axis: i}] for name, a in arrays.items()}, + axis_temperature=axis_temperature, + axis_velocity=axis_velocity, + ) + for i in range(shape[axis]) + ], + axis=axis, + ) + + return na.FunctionArray( + inputs=na.DopplerVectorArray.from_velocity( + velocity=velocity, + wavelength_rest=wavelength_rest, + ), + outputs=result, + ) + + +def _line( + contribution: na.FunctionArray, + temperature: na.AbstractScalar, + mass: na.AbstractScalar, + emission_measure: na.AbstractScalar, + velocity: na.AbstractScalar, + shift: na.AbstractScalar, + width_nonthermal: na.AbstractScalar, + axis_temperature: str, + axis_velocity: str, +) -> na.AbstractScalar: + """ + The spectral radiance of a single line in each bin of velocity. + + The arguments are those of :func:`synthesize`, with any axis along the + lines already indexed away. + """ + width = np.sqrt( + astropy.constants.k_B * temperature / mass + np.square(width_nonthermal) + ) + + # The fraction of the line which falls in each bin, as the difference of + # the cumulative distribution of the profile across the bin. + z = (velocity - shift) / (np.sqrt(2) * width) + cdf = scipy.special.erf(z.to(u.dimensionless_unscaled)) / 2 + fraction = _diff(cdf, axis_velocity) + + radiance = (contribution.outputs * emission_measure * fraction).sum( + axis_temperature + ) + + # The width of each bin in wavelength, which is what turns the radiance in + # a bin into a spectral radiance. + wavelength = (1 + velocity / astropy.constants.c) * contribution.inputs + width_bin = _diff(wavelength, axis_velocity) + + result = radiance / width_bin / (4 * np.pi * u.sr) + + return result.to(u.erg / u.cm**2 / u.s / u.sr / u.AA) + + +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)) + + +def _diff( + a: na.AbstractScalar, + axis: str, +) -> na.AbstractScalar: + """The difference between neighboring vertices along an axis.""" + return a[{axis: slice(1, None)}] - a[{axis: slice(None, -1)}] diff --git a/utu/spectrum/_tests/test_synthesis.py b/utu/spectrum/_tests/test_synthesis.py new file mode 100644 index 0000000..b8389c8 --- /dev/null +++ b/utu/spectrum/_tests/test_synthesis.py @@ -0,0 +1,257 @@ +import astropy.constants +import astropy.units as u +import fiasco +import named_arrays as na +import numpy as np +import pytest + +import utu + +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 _total(result: na.FunctionArray) -> na.AbstractScalar: + """The radiance of each line, summed over the bins of velocity.""" + width = result.inputs.wavelength + width = width[{"velocity": slice(1, None)}] - width[{"velocity": slice(None, -1)}] + return (result.outputs * width).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 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 + + result = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=em, + velocity=velocity, + shift=shift, + ) + + 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 = utu.spectrum.synthesize( + contribution=contribution, + temperature=temperature, + mass=mass, + emission_measure=em[index], + velocity=velocity, + shift=shift[index], + ) + assert np.allclose(pixel.outputs, result.outputs[index], 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) From 89f51a400de2fda879c8698905582a6ba0875a42 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Wed, 14 Oct 2026 20:55:22 -0600 Subject: [PATCH 2/2] Fix the issues the review of `synthesize` found - Raise a ValueError when `velocity` lacks `axis_velocity` or varies with temperature, when `temperature` or `contribution` lacks the axis of temperature, when a temperature is not positive, and when an argument varies along the axis of lines differently than `contribution` does. A velocity axis of the wrong name used to give a result of NaN. - Broadcast arguments with one element along the axis of lines, and return an empty result for a `contribution` with no lines, instead of raising. - Take the fraction of a profile in a bin from the complementary error function of the tail the bin is in. The difference of error functions rounded every bin beyond about 8.5 standard deviations to zero. - Fill a preallocated result one bin of velocity and one block of lines at a time, keeping every intermediate array under `_size_block` elements, and compute the profile once for all lines when nothing it depends on varies along the axis of lines. - Take the width of a bin in wavelength from the width of the bin in velocity, rather than as the difference of two Doppler-shifted wavelengths. - Say in the docs that the emission measure is per bin of temperature, not a DEM per unit log T, and document the n_e n_H convention in `lines()` and `contribution_function()` as well. - Test axis names other than the defaults, every new error, an empty and a broadcast axis of lines, reversed bins, the far wings, nonthermal widths which vary over a map, and synthesis in several blocks. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_017xqpht97bF4FWsACAA4wQa --- docs/tutorials/spectrum.rst | 6 +- utu/spectrum/_lines.py | 17 +- utu/spectrum/_synthesis.py | 266 +++++++++++++++++++------- utu/spectrum/_tests/test_synthesis.py | 172 +++++++++++++++-- 4 files changed, 371 insertions(+), 90 deletions(-) diff --git a/docs/tutorials/spectrum.rst b/docs/tutorials/spectrum.rst index 78dff47..fba3c80 100644 --- a/docs/tutorials/spectrum.rst +++ b/docs/tutorials/spectrum.rst @@ -254,10 +254,6 @@ profiles give the radiance of the line in the table above. .. jupyter-execute:: - wavelength = result.inputs.wavelength - width_bin = ( - wavelength[{"velocity": slice(1, None)}] - - wavelength[{"velocity": slice(None, -1)}] - ) + width_bin = np.diff(result.inputs.wavelength, axis="velocity") (result.outputs * width_bin).sum("velocity") 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 index 26e01e3..5d733f5 100644 --- a/utu/spectrum/_synthesis.py +++ b/utu/spectrum/_synthesis.py @@ -12,6 +12,17 @@ "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, @@ -53,7 +64,9 @@ def synthesize( 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 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 ---------- @@ -67,16 +80,21 @@ def synthesize( :math:`n_e n_H`. temperature The temperatures of the plasma, along `axis_temperature`: the ones - `contribution` was computed at. + `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 at each temperature, - :math:`\int n_e n_H \, dh` through the column in each bin of - temperature, along `axis_temperature`. Every other axis is carried - through to the result, such as the two axes of a map. + 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. @@ -94,16 +112,27 @@ def synthesize( axis_velocity The name of the axis of `velocity`. axis - The name of the axis along the lines of `contribution`. + 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 `emission_measure` but + `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 @@ -113,9 +142,12 @@ def synthesize( left out, which in the transition region and the corona is a small fraction of a milliangstrom, and so is opacity. - The lines are synthesized one at a time, so that the largest array held - at once is the profile of one line at every temperature and in every - pixel, rather than that of every line. + 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 -------- @@ -178,41 +210,132 @@ def synthesize( 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 - arrays = dict( - temperature=_scalar(temperature), - mass=_scalar(mass), - emission_measure=_scalar(emission_measure), - velocity=_scalar(velocity), - shift=_scalar(shift), - width_nonthermal=_scalar(width_nonthermal), + 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), ) - shape = na.shape(contribution.outputs) - - # A single line, such as one picked out of the lines of an ion with an - # integer index, has no axis of lines to loop over. - if axis not in shape: - result = _line( - contribution=contribution, - **arrays, - axis_temperature=axis_temperature, - axis_velocity=axis_velocity, - ) + # 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: - result = na.stack( - [ - _line( - contribution=contribution[{axis: i}], - **{name: a[{axis: i}] for name, a in arrays.items()}, - axis_temperature=axis_temperature, + 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, ) - for i in range(shape[axis]) - ], - axis=axis, - ) + + 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( @@ -223,45 +346,62 @@ def synthesize( ) -def _line( - contribution: na.FunctionArray, +def _fraction( temperature: na.AbstractScalar, mass: na.AbstractScalar, - emission_measure: na.AbstractScalar, velocity: na.AbstractScalar, shift: na.AbstractScalar, width_nonthermal: na.AbstractScalar, - axis_temperature: str, axis_velocity: str, ) -> na.AbstractScalar: """ - The spectral radiance of a single line in each bin of velocity. + 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 arguments are those of :func:`synthesize`, with any axis along the - lines already indexed away. + 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) ) - # The fraction of the line which falls in each bin, as the difference of - # the cumulative distribution of the profile across the bin. z = (velocity - shift) / (np.sqrt(2) * width) - cdf = scipy.special.erf(z.to(u.dimensionless_unscaled)) / 2 - fraction = _diff(cdf, axis_velocity) + z = z.to(u.dimensionless_unscaled).value - radiance = (contribution.outputs * emission_measure * fraction).sum( - axis_temperature - ) + 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) - # The width of each bin in wavelength, which is what turns the radiance in - # a bin into a spectral radiance. - wavelength = (1 + velocity / astropy.constants.c) * contribution.inputs - width_bin = _diff(wavelength, axis_velocity) + erfc0 = scipy.special.erfc(side * z0) + erfc1 = scipy.special.erfc(side * z1) - result = radiance / width_bin / (4 * np.pi * u.sr) + 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. - return result.to(u.erg / u.cm**2 / u.s / u.sr / u.AA) + 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( @@ -269,11 +409,3 @@ def _scalar( ) -> 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)) - - -def _diff( - a: na.AbstractScalar, - axis: str, -) -> na.AbstractScalar: - """The difference between neighboring vertices along an axis.""" - return a[{axis: slice(1, None)}] - a[{axis: slice(None, -1)}] diff --git a/utu/spectrum/_tests/test_synthesis.py b/utu/spectrum/_tests/test_synthesis.py index b8389c8..5746a59 100644 --- a/utu/spectrum/_tests/test_synthesis.py +++ b/utu/spectrum/_tests/test_synthesis.py @@ -4,8 +4,10 @@ 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 @@ -44,11 +46,26 @@ 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.""" - width = result.inputs.wavelength - width = width[{"velocity": slice(1, None)}] - width[{"velocity": slice(None, -1)}] - return (result.outputs * width).sum("velocity") + return _radiance(result).sum("velocity") def _moments(result: na.FunctionArray) -> tuple[na.AbstractScalar, na.AbstractScalar]: @@ -184,8 +201,9 @@ def test_synthesize_single(): def test_synthesize_map(): """ - The axes of the emission measure are carried through, and the shift can - vary along them and with temperature, each pixel on its own. + 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") @@ -193,13 +211,13 @@ def test_synthesize_map(): 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 - result = utu.spectrum.synthesize( - contribution=contribution, - temperature=temperature, - mass=mass, + 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, - velocity=velocity, shift=shift, + width_nonthermal=width, ) assert na.shape(result.outputs) == dict(line=2, velocity=800, x=3, y=2) @@ -207,17 +225,141 @@ def test_synthesize_map(): for i in range(3): for j in range(2): index = {"x": i, "y": j} - pixel = utu.spectrum.synthesize( - contribution=contribution, - temperature=temperature, - mass=mass, + pixel = _synthesize( emission_measure=em[index], - velocity=velocity, 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."""