added a ValueError if trying to extrapolate and update to documentati…

…on (#41)

* added example to illustrate interpolation of edges

* added info about interpolation and custom materials

* added info on energy range for NIST data

* typo on energy range

* change to error if outside energy range

* Update CHANGELOG.rst

* black reformat

* changed ranges to remove now ValueErrors

* Update docs/examples/tutorial1b.rst

Co-authored-by: Shane Maloney <shane.maloney@dias.ie>

Co-authored-by: Shane Maloney <shane.maloney@dias.ie>

CHANGELOG.rst

@@ -1,6 +1,13 @@
 Changelog
 =========
 
+2.1.0 (2022-01-07)
+------------------
+* Improvement to documentation to illustrate custom materials
+* Added tutorial and documentation to illustrate interpolation
+* Going outside of data range (1 keV to 2 MeV) now provides a ValueError instead of failing quietly and filling in zeros for extrapolated values.
+
+
 2.0.1 (2022-10-10)
 ------------------
 * The Material class can now represent substances with multiple and arbitrary constitutients

docs/examples/tutorial1b.rst

@@ -0,0 +1,92 @@
+Edges in the Mass Attenuation Coefficient and Interpolation
+===========================================================
+
+Some elements have complex features in their mass attentuation coefficients.
+The mass attenuation coefficients are interpolated between data points so be careful to include sufficient points to resolve those features if you are looking for high accuracy.
+
+.. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    import astropy.units as u
+    from astropy.visualization import quantity_support
+    quantity_support()
+    from roentgen.absorption import MassAttenuationCoefficient
+
+    cdte_atten = MassAttenuationCoefficient('cdte')
+
+    energy = u.Quantity(np.arange(3, 6, 0.1), 'keV')
+    atten = cdte_atten.func(energy)
+
+    plt.plot(energy, atten)
+    plt.plot(cdte_atten.energy, cdte_atten.data, 'o')
+    plt.yscale('log')
+    plt.xscale('log')
+    plt.ylim(300, 1000)
+    plt.xlim(3, 6)
+    plt.xlabel('Energy [' + str(energy.unit) + ']')
+    plt.ylabel('Mass attenuation Coefficient [' + str(atten.unit) + ']')
+    plt.title(cdte_atten.name + ' undersampling!')
+    plt.show()
+
+The above example is clearly undersampled. Let's add better sampling.
+
+.. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    import astropy.units as u
+    from astropy.visualization import quantity_support
+    quantity_support()
+    from roentgen.absorption import MassAttenuationCoefficient
+
+    cdte_atten = MassAttenuationCoefficient('cdte')
+
+    energy = u.Quantity(np.arange(3, 6, 0.01), 'keV')
+    atten = cdte_atten.func(energy)
+
+    plt.plot(energy, atten)
+    plt.plot(cdte_atten.energy, cdte_atten.data, 'o')
+    plt.yscale('log')
+    plt.xscale('log')
+    plt.ylim(300, 1000)
+    plt.xlim(3, 6)
+    plt.xlabel('Energy [' + str(energy.unit) + ']')
+    plt.ylabel('Mass attenuation Coefficient [' + str(atten.unit) + ']')
+    plt.title(cdte_atten.name + ' better sampling!')
+    plt.show()
+
+This looks much better! Though if we look very closely, we see that we are still undersampling.
+
+.. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    import astropy.units as u
+    from astropy.visualization import quantity_support
+    quantity_support()
+    from roentgen.absorption import MassAttenuationCoefficient
+
+    cdte_atten = MassAttenuationCoefficient('cdte')
+
+    energy = u.Quantity(np.arange(3, 6, 0.01), 'keV')
+    atten = cdte_atten.func(energy)
+
+    plt.plot(energy, atten)
+    plt.plot(cdte_atten.energy, cdte_atten.data, 'o')
+    plt.yscale('log')
+    plt.xscale('log')
+    plt.ylim(600, 900)
+    plt.xlim(3.95, 4.1)
+    plt.xlabel('Energy [' + str(energy.unit) + ']')
+    plt.ylabel('Mass attenuation Coefficient [' + str(atten.unit) + ']')
+    plt.title(cdte_atten.name + ' still undersampled')
+    plt.show()
+
+For many calculations, this small difference may not matter.

docs/guide/transmission.rst

@@ -5,13 +5,14 @@ The purpose of this guide is to present an overview of the `roentgen.absorption`
 Mass Attenuation Coefficient
 ----------------------------
 The primary component that mediates the x-ray attenuation through a material is its mass attenuation coefficient.
-These tabulated values can be inspected using `roentgen.absorption.MassAttenuationCoefficient`.
+These tabulated values can be inspected using the `roentgen.absorption.MassAttenuationCoefficient` object.
 To create one::
 
     from roentgen.absorption import MassAttenuationCoefficient
     si_matten = MassAttenuationCoefficient('Si')
 
 Tabulated values for all elements are provided as well as additional specialized materials.
+These data are provided by the U.S National Institute of Standards and Technology (`NIST <https://physics.nist.gov/PhysRefData/XrayMassCoef/tab3.html>`__) and range from 1 keV to 20 MeV.
 Elements can be specificied by their symbol name or by their full name (e.g. Si, Silicon).
 A list of all of the elements is provided by::
 
@@ -21,7 +22,7 @@ Specialized materials, referred to as compounds, are also available. A complete
 
     roentgen.compounds
 
-Here is the mass attenuation coefficient for Silicon.
+Here is the mass attenuation coefficient data for Silicon.
 
 .. plot::
     :include-source:
@@ -31,13 +32,41 @@ Here is the mass attenuation coefficient for Silicon.
     from roentgen.absorption import MassAttenuationCoefficient
 
     si_matten = MassAttenuationCoefficient('Si')
-    plt.plot(si_matten.energy, si_matten.data)
+    plt.plot(si_matten.energy, si_matten.data, 'o-')
     plt.yscale('log')
     plt.xscale('log')
     plt.xlabel('Energy [' + str(si_matten.energy[0].unit) + ']')
     plt.ylabel('Mass Attenuation Coefficient [' + str(si_matten.data[0].unit) + ']')
     plt.title(si_matten.name)
 
+You can get the interpolated value at any other point
+
+.. plot::
+    :include-source:
+
+    import astropy.units as u
+    import matplotlib.pyplot as plt
+    from roentgen.absorption import MassAttenuationCoefficient
+
+    si_matten = MassAttenuationCoefficient('Si')
+
+    energy = u.Quantity(np.arange(1, 10, 0.01), 'keV')
+    interpol_atten = si_matten.func(energy)
+
+    plt.plot(energy, interpol_atten, 'x-', label='Interpolated points')
+    plt.plot(si_matten.energy, si_matten.data, 'o', label='Data')
+
+    plt.xlim(1, 10)
+    plt.ylim(10, 5000)
+    plt.yscale('log')
+    plt.xscale('log')
+    plt.xlabel('Energy [' + str(si_matten.energy[0].unit) + ']')
+    plt.ylabel('Mass Attenuation Coefficient [' + str(si_matten.data[0].unit) + ']')
+    plt.title(si_matten.name)
+    plt.legend()
+
+Be careful that you interpolate with sufficient resolution to make out complex edges if high accuracy is required.
+
 Material
 --------
 In order to determine the x-ray attenuation through a material the `roentgen.absorption.Material` object is provided.
@@ -128,6 +157,21 @@ As a simple example, here is the transmission of x-rays through 10 meters of air
 This plot shows that air, though not a dense material, does block low energy x-rays over long distances.
 For convenience, the function `~roentgen.util.density_ideal_gas` is provided which can calculate the density of a gas given a pressure and temperature.
 
+It is also possible to create custom materials using a combination of elements and compounds::
+
+    >>> from roentgen.absorption import Material
+    >>> steel = Material({"Fe": 0.98, "C": 0.02}, 1 * u.cm)
+    >>> bronze = Material({"Cu": 0.88, "Sn": 0.12}, 1 * u.cm)
+    >>> salt_water = Material({"water": 0.97, "Na": 0.015, "Cl": 0.015}, 1 * u.cm)
+
+The fractions need not be normalized. It will normalize them for you.
+The density will be calculated automatically using the known densities but this is likely not a good assumption so you should provide your own density::
+    
+    >>> bronze = Material({"Cu": 0.88, "Sn": 0.12}, 1 * u.cm)
+    >>> bronze.density
+    <Quantity 7746.6 kg / m3>
+    >>> bronze = Material({"Cu": 0.88, "Sn": 0.12}, 1 * u.cm, density=8.73 * u.g u.cm**-3)
+
 Stack
 -----
 Materials can be added together to form more complex optical paths.

roentgen/absorption/material.py

@@ -160,6 +160,11 @@ def transmission(self, energy):
         ----------
         energy : `astropy.units.Quantity`
             An array of energies in keV
+
+        Raises
+        ------
+        ValueError
+            If energy is outside of the interpolation range of 1 keV to 20 MeV.
         """
         coefficients = self.mass_attenuation_coefficient(energy)
         transmission = np.exp(-coefficients * self.density * self.thickness)
@@ -173,6 +178,11 @@ def absorption(self, energy):
         ----------
         energy : `astropy.units.Quantity`
             An array of energies in keV.
+
+        Raises
+        ------
+        ValueError
+            If energy is outside of the interpolation range of 1 keV to 20 MeV.
         """
         return 1.0 - self.transmission(energy)
 
@@ -185,6 +195,11 @@ def linear_attenuation_coefficient(self, energy: u.keV):
         ----------
         energy : `astropy.units.Quantity`
             An array of energies in keV.
+
+        Raises
+        ------
+        ValueError
+            If energy is outside of the interpolation range of 1 keV to 20 MeV.
         """
         return self.mass_attenuation_coefficient(energy) * self.density
 
@@ -233,6 +248,11 @@ def transmission(self, energy):
         ----------
         energy : `astropy.units.Quantity`
             An array of energies in keV
+
+        Raises
+        ------
+        ValueError
+            If energy is outside of the interpolation range of 1 keV to 20 MeV.
         """
         transmission = np.ones(len(energy), dtype=float)
         for material in self.materials:
@@ -247,6 +267,11 @@ def absorption(self, energy):
         ----------
         energy : `astropy.units.Quantity`
             An array of energies in keV.
+
+        Raises
+        ------
+        ValueError
+            If energy is outside of the interpolation range of 1 keV to 20 MeV.
         """
         return 1.0 - self.transmission(energy)
 
@@ -313,6 +338,11 @@ def response(self, energy):
         ----------
         energy : `astropy.units.Quantity`
             An array of energies in keV.
+
+        Raises
+        ------
+        ValueError
+            If energy is outside of the interpolation range of 1 keV to 20 MeV.
         """
         # calculate the transmission
         transmission = np.ones(len(energy), dtype=float)
@@ -348,6 +378,8 @@ class MassAttenuationCoefficient(object):
     func : `lambda func`
         A function which returns the interpolated mass attenuation value at
         any given energy. Energies must be given by an `astropy.units.Quantity`.
+        The interpolation range is 1 keV to 20 MeV.
+        Going outside that range will result in a ValueError.
 
     """
 
@@ -390,8 +422,7 @@ def __init__(self, material):
         self._f = interpolate.interp1d(
             data_energy_kev,
             data_attenuation_coeff,
-            bounds_error=False,
-            fill_value=0.0,
+            bounds_error=True,
             assume_sorted=True,
         )
         self.func = lambda x: u.Quantity(10 ** self._f(np.log10(x.to("keV").value)), "cm^2/g")

roentgen/data/README.rst

@@ -21,6 +21,7 @@ elements.csv
 ------------
 Source: `https://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html <https://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html>`__
 Provides translation between atomic number, element names, and element symbols.
+Energy ranges from 1 keV to 20 MeV.
 
 - **z**: Atomic number
 - **name**: Element name
@@ -33,6 +34,7 @@ compounds_mixtures.csv
 Source: `https://physics.nist.gov/PhysRefData/XrayMassCoef/tab2.html <https://physics.nist.gov/PhysRefData/XrayMassCoef/tab2.html>`__
 Provides a list of all compounds and mixtures supported
 density is in units of g/cm^3
+Energy ranges from 1 keV to 20 MeV.
 
 emission_lines.csv
 ------------------

roentgen/tests/test_material.py

@@ -50,11 +50,11 @@ def test_mass_atten(mass_atten):
 
 
 def test_returns_quantity(mass_atten):
-    assert isinstance(mass_atten.func(1 * u.keV), u.Quantity)
+    assert isinstance(mass_atten.func(2 * u.keV), u.Quantity)
 
 
 def test_number_of_energies(mass_atten):
-    energy = u.Quantity(np.arange(1, 1000), "keV")
+    energy = u.Quantity(np.arange(2, 1000), "keV")
     atten = mass_atten.func(energy)
     assert len(energy) == len(atten)
 
@@ -95,7 +95,7 @@ def thick_material(request):
 
 def test_opaque(thick_material):
     # check that extremely large amounts of material mean no transmission
-    assert thick_material.transmission(1 * u.keV) < 1e-6
+    assert thick_material.transmission(2 * u.keV) < 1e-6
 
 
 @pytest.fixture(params=all_materials)
@@ -105,7 +105,7 @@ def thin_material(request):
 
 def test_transparent(thin_material):
     # check that extremely large amounts of material mean no transmission
-    assert thin_material.transmission(1 * u.keV) > 0.90
+    assert thin_material.transmission(2 * u.keV) > 0.90
 
 
 @pytest.mark.parametrize(
@@ -219,3 +219,23 @@ def test_density_calculation():
     this_mat = Material(a, 5 * u.m)
     calc_density = get_material_density("Fe") * 0.5 + get_material_density("C") * 0.5
     assert np.isclose(this_mat.density, calc_density)
+
+
+def test_raise_outside_of_data_range():
+    """Test that ValueError is raised is trying to get values outside of data range 1 keV to 20 MeV."""
+    mat = Material("Fe", 1 * u.m)
+    energy = u.Quantity(np.arange(0.1, 10, 0.1), "keV")
+    # below 1 keV
+    with pytest.raises(ValueError):
+        mat.absorption(energy)
+
+    with pytest.raises(ValueError):
+        mat.transmission(energy)
+
+    # above 20 MeV
+    energy = u.Quantity(np.arange(10, 23, 0.1), "MeV")
+    with pytest.raises(ValueError):
+        mat.absorption(energy)
+
+    with pytest.raises(ValueError):
+        mat.transmission(energy)

roentgen/tests/test_response.py

@@ -42,3 +42,18 @@ def test_optical_path_stack_input():
 def test_bad_optical_path():
     with pytest.raises(TypeError):
         Response(optical_path="Si", detector=thin_material)
+
+
+def test_raise_outside_of_data_range():
+    """Test that ValueError is raised is trying to get values outside of data range 1 keV to 20 MeV."""
+    resp = Response(optical_path=thin_material, detector=thin_material)
+
+    energy = u.Quantity(np.arange(0.1, 10, 0.1), "keV")
+    # below 1 keV
+    with pytest.raises(ValueError):
+        resp.response(energy)
+
+    # above 20 MeV
+    energy = u.Quantity(np.arange(10, 23, 0.1), "MeV")
+    with pytest.raises(ValueError):
+        resp.response(energy)