From d9c5953fa2cc9d2bac7e3b116b75e59d3352d9d6 Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Wed, 23 Oct 2024 17:05:25 -0400
Subject: [PATCH] Use log-log interpolation for xray f''. Fixes #55.

---
 doc/sphinx/guide/using.rst |  2 +-
 periodictable/xsf.py       | 31 ++++++++++++++++---------------
 2 files changed, 17 insertions(+), 16 deletions(-)

diff --git a/doc/sphinx/guide/using.rst b/doc/sphinx/guide/using.rst
index ad45a35..7c5bf5d 100644
--- a/doc/sphinx/guide/using.rst
+++ b/doc/sphinx/guide/using.rst
@@ -29,7 +29,7 @@ Access particular elements as symbols:
     B absorption 767.0
     >>> print("Ni f1/f2 for Cu K-alpha X-rays f'=%.5f f''=%.5f"
     ...       % Ni.xray.scattering_factors(wavelength=Cu.K_alpha))
-    Ni f1/f2 for Cu K-alpha X-rays f'=25.02293 f''=0.52493
+    Ni f1/f2 for Cu K-alpha X-rays f'=25.02293 f''=0.52491
 
 Access isotopes using mass number subscripts:
 
diff --git a/periodictable/xsf.py b/periodictable/xsf.py
index 7a77d63..eb52b92 100644
--- a/periodictable/xsf.py
+++ b/periodictable/xsf.py
@@ -178,7 +178,6 @@
 .. [#Deslattes2003] R. D. Deslattes, E. G. Kessler, Jr., P. Indelicato, L. de Billy,
        E. Lindroth, and J. Anton.  Rev. Mod. Phys. 75, 35-99 (2003).
 """
-from __future__ import with_statement
 __all__ = ['Xray', 'init', 'init_spectral_lines',
            'xray_energy', 'xray_wavelength',
            'xray_sld', 'xray_sld_from_atoms',
@@ -187,7 +186,7 @@
            ]
 import os.path
 
-import numpy
+import numpy as np
 from numpy import nan, pi, exp, sin, cos, sqrt, radians
 
 from .core import Element, Ion, default_table, get_data_path
@@ -218,7 +217,7 @@ def xray_wavelength(energy):
 
         $c$ = speed of light in m/s
     """
-    return planck_constant/electron_volt*speed_of_light/numpy.asarray(energy)*1e7
+    return planck_constant/electron_volt*speed_of_light/np.asarray(energy)*1e7
 
 def xray_energy(wavelength):
     r"""
@@ -242,7 +241,7 @@ def xray_energy(wavelength):
 
         $c$ = speed of light in m/s
     """
-    return planck_constant/electron_volt*speed_of_light/numpy.asarray(wavelength)*1e7
+    return planck_constant/electron_volt*speed_of_light/np.asarray(wavelength)*1e7
 
 class Xray(object):
     """
@@ -265,8 +264,8 @@ def _gettable(self):
             filename = os.path.join(self._nff_path,
                                     self.element.symbol.lower()+".nff")
             if self.element.symbol != 'n' and os.path.exists(filename):
-                xsf = numpy.loadtxt(filename, skiprows=1).T
-                xsf[1, xsf[1] == -9999.] = numpy.nan
+                xsf = np.loadtxt(filename, skiprows=1).T
+                xsf[1, xsf[1] == -9999.] = np.nan
                 xsf[0] *= 0.001  # Use keV in table rather than eV
                 self._table = xsf
         return self._table
@@ -298,11 +297,13 @@ def scattering_factors(self, energy=None, wavelength=None):
         if energy is None:
             raise TypeError('X-ray scattering factors need wavelength or energy')
 
-        scalar = numpy.isscalar(energy)
+        scalar = np.isscalar(energy)
         if scalar:
-            energy = numpy.array([energy])
-        f1 = numpy.interp(energy, xsf[0], xsf[1], left=nan, right=nan)
-        f2 = numpy.interp(energy, xsf[0], xsf[2], left=nan, right=nan)
+            energy = np.array([energy])
+        f1 = np.interp(energy, xsf[0], xsf[1], left=nan, right=nan)
+        f2 = np.exp(np.interp(
+            np.log(energy), np.log(xsf[0]), np.log(xsf[2]),
+            left=nan, right=nan))
         if scalar:
             f1, f2 = f1[0], f2[0]
         return f1, f2
@@ -505,10 +506,10 @@ def mirror_reflectivity(compound, density=None, natural_density=None,
         wavelength = xray_wavelength(energy)
     assert wavelength is not None, "scattering calculation needs energy or wavelength"
     angle = radians(angle)
-    if numpy.isscalar(wavelength):
-        wavelength = numpy.array([wavelength])
-    if numpy.isscalar(angle):
-        angle = numpy.array([angle])
+    if np.isscalar(wavelength):
+        wavelength = np.array([wavelength])
+    if np.isscalar(angle):
+        angle = np.array([angle])
     nv = index_of_refraction(compound=compound,
                              density=density, natural_density=natural_density,
                              wavelength=wavelength)
@@ -701,7 +702,7 @@ def sld_table(wavelength=None, table=None):
          He   1.03   0.00
          Li   3.92   0.00
          Be  13.93   0.01
-          B  18.40   0.02
+          B  18.40   0.01
           C  18.71   0.03
           N   6.88   0.02
           O   9.74   0.04