Merge pull request #25 from ehsteve/bug_fix

Bug fix

docs/examples/tutorial1.rst

@@ -4,11 +4,14 @@ Getting and plotting the Mass Attenuation Coefficient
 .. plot::
     :include-source:
 
-    from roentgen.absorption import MassAttenuationCoefficient
     import numpy as np
-    import astropy.units as u
     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(1, 1000), 'keV')

docs/examples/tutorial2.rst

@@ -8,6 +8,9 @@ Finding the x-ray absorption of different detector materials
     import matplotlib.pyplot as plt
 
     import astropy.units as u
+    from astropy.visualization import quantity_support
+    quantity_support()
+
     from roentgen.absorption import Material
 
     thickness = 500 * u.micron

docs/examples/tutorial5.rst

@@ -4,11 +4,15 @@ Finding the response of an x-ray detector
 .. plot::
     :include-source:
 
-    from roentgen.absorption import Material, Response
     import numpy as np
-    import astropy.units as u
     import matplotlib.pyplot as plt
 
+    import astropy.units as u
+    from astropy.visualization import quantity_support
+    quantity_support()
+
+    from roentgen.absorption import Material, Response
+
     optical_path = [Material('air', 2 * u.m), Material('mylar', 5 * u.micron), Material('Al', 5 * u.micron)]
     detector = Material('Si', 500 * u.micron)
 

docs/examples/tutorial7.rst

@@ -0,0 +1,26 @@
+Checking on edges in the Mass Attenuation Coefficient
+=====================================================
+
+.. 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
+
+    e = np.linspace(1e-3, 1e-2, 1000)*u.MeV
+    te = MassAttenuationCoefficient('Te')
+    plt.plot(e, te.func(e), 'x', label='Interpolated')
+    plt.plot(te.energy, te.data, color='black', label='Data')
+    plt.loglog()
+    plt.xlim(4e-3, 6e-3)
+    plt.ylim(1e2, 1e3)
+    plt.xlabel('Photon Energy (MeV)')
+    plt.ylabel(r'$\mu / \rho$ (cm$^{2}$ g$^{-1}$)')
+    plt.title('Te')
+    plt.legend()
+    plt.show()

roentgen/absorption/material.py

@@ -309,12 +309,25 @@ def __init__(self, material):
         self.energy = u.Quantity(data[:, 0] * 1000, "keV")
         self.data = u.Quantity(data[:, 1], "cm^2/g")
 
+        self._remove_double_vals_from_data()
+
         data_energy_kev = np.log10(self.energy.value)
         data_attenuation_coeff = np.log10(self.data.value)
         self._f = interpolate.interp1d(
             data_energy_kev, data_attenuation_coeff, bounds_error=False,
-            fill_value=0.0
+            fill_value=0.0, assume_sorted=True
         )
         self.func = lambda x: u.Quantity(
             10 ** self._f(np.log10(x.to("keV").value)), "cm^2/g"
         )
+
+    def _remove_double_vals_from_data(self):
+        """Remove double-values energy values. Edges are represented with
+        the same energy index and at the bottom and top value of the edge. This
+        must be removed to enable correct interpolation."""
+        uniq, count = np.unique(self.energy, return_counts=True)
+        duplicates = uniq[count > 1]
+        for this_dup in duplicates:
+            ind = (self.energy == this_dup).nonzero()
+            # shift the first instance of the energy, the bottom of the edge
+            self.energy[ind[0][0]] -= 1e-3 * u.eV

roentgen/tests/test_massatten.py

@@ -0,0 +1,18 @@
+import pytest
+
+from numpy.testing import assert_allclose
+
+import astropy.units as u
+
+from roentgen.absorption.material import MassAttenuationCoefficient
+
+
+def test_interpolate_matches_at_data():
+    for this_element in ['Te', 'Si', 'Ge', 'cdte']:
+        te = MassAttenuationCoefficient(this_element)
+        assert_allclose(te.data, te.func(te.energy))
+
+        # now change one point and make sure it no longer works
+        te.energy[0] = 10 * u.keV
+        with pytest.raises(AssertionError):
+            assert_allclose(te.data, te.func(te.energy))