Expose objcryst's anisotropic Pseudo-Voigt reflection profile

news/expose-aniso-pv.rst

@@ -0,0 +1,25 @@
+**Added:**
+
+* Add explicitly constructable reflection profiles
+* Expose objcryst's anisotropic Pseudo-Voigt profile
+* Add PowderPatternDiffraction.SetProfile
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/extensions/powderpatterndiffraction_ext.cpp

@@ -25,12 +25,24 @@
 
 #include <ObjCryst/ObjCryst/General.h>
 #include <ObjCryst/ObjCryst/PowderPattern.h>
+#include <ObjCryst/ObjCryst/ReflectionProfile.h>
 #include <ObjCryst/ObjCryst/ScatteringData.h>
 
 namespace bp = boost::python;
 using namespace boost::python;
 using namespace ObjCryst;
 
+namespace
+{
+
+void SetProfileCopy(PowderPatternDiffraction &pdiff,
+                    const ReflectionProfile &profile)
+{
+    pdiff.SetProfile(profile.CreateCopy());
+}
+
+} // namespace
+
 
 void wrap_powderpatterndiffraction()
 {
@@ -60,6 +72,8 @@ void wrap_powderpatterndiffraction()
                 (ReflectionProfile& (PowderPatternDiffraction::*)())
                 &PowderPatternDiffraction::GetProfile,
                 return_internal_reference<>())
+        .def("SetProfile", &SetProfileCopy, bp::arg("profile"),
+                "Install an independent copy of a reflection profile.")
         .def("SetExtractionMode",
                 &PowderPatternDiffraction::SetExtractionMode,
                 (bp::arg("extract")=true, bp::arg("init")=false))

src/extensions/reflectionprofile_ext.cpp

@@ -17,7 +17,10 @@
  *
  *****************************************************************************/
 
+#define BOOST_PYTHON_MAX_ARITY 20
+
 #include <boost/python/class.hpp>
+#include <boost/python/init.hpp>
 #include <boost/python/manage_new_object.hpp>
 #include <boost/python/pure_virtual.hpp>
 #undef B0
@@ -113,6 +116,8 @@ void wrap_reflectionprofile()
              (bp::arg("relativeIntensity"), bp::arg("xcenter"),
               bp::arg("h"), bp::arg("k"), bp::arg("l")),
               "Return the full profile width at a given relative intensity for reflection (h, k, l) around `xcenter`.")
+        .def("IsAnisotropic", &ReflectionProfile::IsAnisotropic,
+             "Return whether the profile depends on the reflection indices.")
         .def("XMLOutput",
              pure_virtual((void (ReflectionProfile::*)(ostream &, int) const) & ReflectionProfile::XMLOutput),
              (bp::arg("os"), bp::arg("indent")),
@@ -121,4 +126,44 @@ void wrap_reflectionprofile()
              pure_virtual((void (ReflectionProfile::*)(istream &, const XMLCrystTag &))&ReflectionProfile::XMLInput),
              (bp::arg("is"), bp::arg("tag")),
              "Load ReflectionProfile parameters from an XML stream and tag.");
+
+    class_<ReflectionProfilePseudoVoigt, bases<ReflectionProfile> >(
+        "ReflectionProfilePseudoVoigt", init<>())
+        .def(init<const ReflectionProfilePseudoVoigt &>(bp::arg("old")))
+        .def("CreateCopy", &ReflectionProfilePseudoVoigt::CreateCopy,
+             return_value_policy<manage_new_object>(),
+             "Return an independent copy of this profile.")
+        .def("SetProfilePar", &ReflectionProfilePseudoVoigt::SetProfilePar,
+             (bp::arg("fwhmCagliotiW"), bp::arg("fwhmCagliotiU") = 0,
+              bp::arg("fwhmCagliotiV") = 0, bp::arg("eta0") = 0.5,
+              bp::arg("eta1") = 0),
+             "Set the isotropic pseudo-Voigt profile parameters.");
+
+    class_<ReflectionProfilePseudoVoigtAnisotropic,
+           bases<ReflectionProfile> >(
+        "ReflectionProfilePseudoVoigtAnisotropic", init<>())
+        .def(init<const ReflectionProfilePseudoVoigtAnisotropic &>(
+            bp::arg("old")))
+        .def("CreateCopy",
+             &ReflectionProfilePseudoVoigtAnisotropic::CreateCopy,
+             return_value_policy<manage_new_object>(),
+             "Return an independent copy of this profile.")
+        .def(
+            "SetProfilePar",
+            &ReflectionProfilePseudoVoigtAnisotropic::SetProfilePar,
+            (bp::arg("fwhmCagliotiW"), bp::arg("fwhmCagliotiU") = 0,
+             bp::arg("fwhmCagliotiV") = 0, bp::arg("fwhmGaussP") = 0,
+             bp::arg("fwhmLorentzX") = 0,
+             bp::arg("fwhmLorentzY") = 0,
+             bp::arg("fwhmLorentzGammaHH") = 0,
+             bp::arg("fwhmLorentzGammaKK") = 0,
+             bp::arg("fwhmLorentzGammaLL") = 0,
+             bp::arg("fwhmLorentzGammaHK") = 0,
+             bp::arg("fwhmLorentzGammaHL") = 0,
+             bp::arg("fwhmLorentzGammaKL") = 0,
+             bp::arg("pseudoVoigtEta0") = 0,
+             bp::arg("pseudoVoigtEta1") = 0,
+             bp::arg("asymA0") = 1, bp::arg("asymA1") = 0,
+             bp::arg("asymA2") = 0),
+            "Set the anisotropic pseudo-Voigt profile parameters.");
 }

src/pyobjcryst/reflectionprofile.py

@@ -20,6 +20,22 @@
 ``XMLOutput`` / ``XMLInput`` and ``CreateCopy``. ``GetProfile``
 accepts a Python sequence or numpy array for ``x``.
 
+Concrete isotropic and anisotropic pseudo-Voigt profiles can be installed
+with ``PowderPatternDiffraction.SetProfile``. The diffraction component
+stores an independent copy, so one configured profile can be reused as a
+template for multiple phases::
+
+    from pyobjcryst.reflectionprofile import (
+        ReflectionProfilePseudoVoigtAnisotropic,
+    )
+
+    profile = ReflectionProfilePseudoVoigtAnisotropic()
+    profile.GetPar("W").SetValue(1e-6)
+    profile.GetPar("G_HH").SetValue(2e-6)
+
+    for pdiff in powder_pattern.get_crystalline_components():
+        pdiff.SetProfile(profile)
+
 Example
 -------
 
@@ -48,6 +64,16 @@
     y_broader = profile.GetProfile(window, xcenter, 1, 0, 0)
 """
 
-__all__ = ["ReflectionProfile", "ReflectionProfileType"]
+__all__ = [
+    "ReflectionProfile",
+    "ReflectionProfilePseudoVoigt",
+    "ReflectionProfilePseudoVoigtAnisotropic",
+    "ReflectionProfileType",
+]
 
-from pyobjcryst._pyobjcryst import ReflectionProfile, ReflectionProfileType
+from pyobjcryst._pyobjcryst import (
+    ReflectionProfile,
+    ReflectionProfilePseudoVoigt,
+    ReflectionProfilePseudoVoigtAnisotropic,
+    ReflectionProfileType,
+)

tests/test_reflectionprofile.py

@@ -7,6 +7,11 @@
 
 from pyobjcryst.powderpattern import PowderPattern
 from pyobjcryst.refinableobj import RefinableObj
+from pyobjcryst.reflectionprofile import (
+    ReflectionProfile,
+    ReflectionProfilePseudoVoigt,
+    ReflectionProfilePseudoVoigtAnisotropic,
+)
 
 
 class TestReflectionProfile(unittest.TestCase):
@@ -19,17 +24,120 @@ def prepare_fixture(self, loadcifdata):
     def setUp(self):
         """Set up a ReflectionProfile instance for testing."""
         x = np.linspace(0, 40, 1000)
-        c = self.loadcifdata("paracetamol.cif")
+        self.crystal = self.loadcifdata("paracetamol.cif")
 
         self.pp = PowderPattern()
         self.pp.SetWavelength(0.7)
         self.pp.SetPowderPatternX(np.deg2rad(x))
         self.pp.SetPowderPatternObs(np.ones_like(x))
 
-        self.ppd = self.pp.AddPowderPatternDiffraction(c)
+        self.ppd = self.pp.AddPowderPatternDiffraction(self.crystal)
 
         self.profile = self.ppd.GetProfile()
 
+    def test_concrete_pseudo_voigt_profiles(self):
+        """Concrete isotropic and anisotropic profiles are
+        constructible."""
+        isotropic = ReflectionProfilePseudoVoigt()
+        anisotropic = ReflectionProfilePseudoVoigtAnisotropic()
+
+        self.assertIsInstance(isotropic, ReflectionProfile)
+        self.assertIsInstance(anisotropic, ReflectionProfile)
+        self.assertFalse(isotropic.IsAnisotropic())
+        self.assertTrue(anisotropic.IsAnisotropic())
+
+        expected = {
+            "U",
+            "V",
+            "W",
+            "P",
+            "X",
+            "Y",
+            "G_HH",
+            "G_KK",
+            "G_LL",
+            "G_HK",
+            "G_HL",
+            "G_KL",
+            "Eta0",
+            "Eta1",
+            "Asym0",
+            "Asym1",
+            "Asym2",
+        }
+        self.assertEqual(
+            expected,
+            {
+                anisotropic.GetPar(i).GetName()
+                for i in range(anisotropic.GetNbPar())
+            },
+        )
+
+    def test_anisotropic_set_profile_par(self):
+        """SetProfilePar assigns widths and retains symmetric default
+        asymmetry."""
+        profile = ReflectionProfilePseudoVoigtAnisotropic()
+        profile.SetProfilePar(1e-6, fwhmLorentzX=2e-6)
+        profile.GetPar("P").SetValue(3e-6)
+
+        self.assertAlmostEqual(profile.GetPar("W").GetValue(), 1e-6)
+        self.assertAlmostEqual(profile.GetPar("X").GetValue(), 2e-6)
+        self.assertAlmostEqual(profile.GetPar("P").GetValue(), 3e-6)
+        self.assertAlmostEqual(profile.GetPar("Asym0").GetValue(), 1.0)
+
+    def test_anisotropic_profile_depends_on_hkl(self):
+        """Anisotropic Lorentz coefficients produce direction
+        dependence."""
+        profile = ReflectionProfilePseudoVoigtAnisotropic()
+        profile.SetProfilePar(
+            1e-6,
+            fwhmLorentzGammaHH=1e-3,
+            fwhmLorentzGammaKK=2e-3,
+            pseudoVoigtEta0=1,
+        )
+        center = 0.5
+        x = np.linspace(center - 0.1, center + 0.1, 401)
+
+        h00 = profile.GetProfile(x, center, 1, 0, 0)
+        zero_k0 = profile.GetProfile(x, center, 0, 1, 0)
+
+        self.assertFalse(np.allclose(h00, zero_k0))
+
+    def test_set_profile_copies_reusable_template(self):
+        """One profile template can safely initialize multiple
+        phases."""
+        second = self.pp.AddPowderPatternDiffraction(self.crystal)
+        template = ReflectionProfilePseudoVoigtAnisotropic()
+        template.GetPar("W").SetValue(1e-6)
+        template.GetPar("G_HH").SetValue(2e-6)
+
+        for pdiff in (self.ppd, second):
+            pdiff.SetProfile(template)
+
+        first_profile = self.ppd.GetProfile()
+        second_profile = second.GetProfile()
+        self.assertIsInstance(
+            first_profile, ReflectionProfilePseudoVoigtAnisotropic
+        )
+        self.assertIsInstance(
+            second_profile, ReflectionProfilePseudoVoigtAnisotropic
+        )
+        self.assertAlmostEqual(first_profile.GetPar("W").GetValue(), 1e-6)
+        self.assertAlmostEqual(second_profile.GetPar("G_HH").GetValue(), 2e-6)
+
+        first_profile.GetPar("W").SetValue(3e-6)
+        self.assertAlmostEqual(template.GetPar("W").GetValue(), 1e-6)
+        self.assertAlmostEqual(second_profile.GetPar("W").GetValue(), 1e-6)
+
+    def test_set_profile_accepts_temporary(self):
+        """A temporary concrete profile is safely copied into the
+        phase."""
+        self.ppd.SetProfile(ReflectionProfilePseudoVoigt())
+        profile = self.ppd.GetProfile()
+
+        self.assertIsInstance(profile, ReflectionProfilePseudoVoigt)
+        self.assertFalse(profile.IsAnisotropic())
+
     def test_get_computed_profile(self):
         """Sample a profile slice and verify broadening lowers the peak
         height."""