diff --git a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp
@@ -86,7 +86,9 @@ fetchSymmetryOperations(const ObjCryst::SpaceGroup& spacegroup)
     assert(nbtran * last <= nbsym);
     for (int nt = 0; nt < nbtran; ++nt)
     {
-        const REAL* pt = sgtrans[nt].tr;
+        // Keep this compatible with ObjCryst builds where REAL is either
+        // float or double.
+        const auto* pt = sgtrans[nt].tr;
         R3::Vector sgt(pt[0], pt[1], pt[2]);
         for (int i = 0; i < last; ++i)
         {

diff --git a/src/diffpy/srreal/PDFCalculator.cpp b/src/diffpy/srreal/PDFCalculator.cpp
@@ -187,10 +187,7 @@ QuantityType PDFCalculator::getExtendedPDF() const
 QuantityType PDFCalculator::getExtendedRDF() const
 {
     QuantityType rdf(this->countExtendedPoints());
-    const double& totocc = mstructure_cache.totaloccupancy;
-    double sfavg = this->sfAverage();
-    double rdf_scale = (totocc * sfavg == 0.0) ? 0.0 :
-        1.0 / (totocc * sfavg * sfavg);
+    double rdf_scale = this->getRDFScale();
     QuantityType::iterator iirdf = rdf.begin();
     QuantityType::const_iterator iival, iival_last;
     iival = this->value().begin() +
@@ -491,9 +488,7 @@ void PDFCalculator::resetValue()
     // when applicable, configure linear baseline
     if (this->getBaseline()->type() == "linear")
     {
-        double partialpdfscale =
-            (0.0 == mstructure_cache.totaloccupancy) ? 0.0 :
-            mstructure_cache.activeoccupancy / mstructure_cache.totaloccupancy;
+        double partialpdfscale = this->getPartialPDFScale();
         double pnumdensity = partialpdfscale * mstructure->numberDensity();
         PDFBaseline& bl = *(this->getBaseline());
         bl.setDoubleAttr("slope", -4 * M_PI * pnumdensity);
@@ -677,6 +672,24 @@ double PDFCalculator::sfAverage() const
 }
 
 
+double PDFCalculator::getPartialPDFScale() const
+{
+    const double totocc = mstructure_cache.totaloccupancy;
+    return (totocc == 0.0) ? 0.0 :
+        (mstructure_cache.activeoccupancy / totocc);
+}
+
+
+double PDFCalculator::getRDFScale() const
+{
+    const double& totocc = mstructure_cache.totaloccupancy;
+    double sfavg = this->sfAverage();
+    double rv = (totocc * sfavg == 0.0) ? 0.0 :
+        1.0 / (totocc * sfavg * sfavg);
+    return rv;
+}
+
+
 void PDFCalculator::cacheStructureData()
 {
     int cntsites = this->countSites();

diff --git a/src/diffpy/srreal/PDFCalculator.hpp b/src/diffpy/srreal/PDFCalculator.hpp
@@ -114,6 +114,14 @@ class PDFCalculator :
         // support for PQEvaluatorOptimized
         virtual void stashPartialValue();
         virtual void restorePartialValue();
+        /// activeoccupancy / totaloccupancy used by baseline background term
+        double getPartialPDFScale() const;
+        /// RDF scale from the cached total occupancy and average scattering factor
+        double getRDFScale() const;
+        /// effective scattering factor at a given site scaled by occupancy
+        const double& sfSite(int) const;
+        /// average scattering factor
+        double sfAverage() const;
 
     private:
 
@@ -145,10 +153,6 @@ class PDFCalculator :
         void cutRipplePoints(QuantityType& y) const;
 
         // structure factors - fast lookup by site index
-        /// effective scattering factor at a given site scaled by occupancy
-        const double& sfSite(int) const;
-        /// average scattering factor
-        double sfAverage() const;
         void cacheStructureData();
         void cacheRlimitsData();
 

diff --git a/src/diffpy/srreal/R3linalg.cpp b/src/diffpy/srreal/R3linalg.cpp
@@ -16,6 +16,8 @@
 *
 *****************************************************************************/
 
+#include <algorithm>
+#include <cmath>
 #include <boost/functional/hash.hpp>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
@@ -82,6 +84,88 @@ const Matrix& inverse(const Matrix& A)
 }
 
 
+void eigen_solve_3x3(const Matrix& A, Vector& w, Matrix& V)
+{
+    V = identity();
+    Matrix m = A;
+
+    const int max_iter = 50;
+    const double eps = 1e-10;
+
+    for (int iter = 0; iter < max_iter; ++iter)
+    {
+        double max_off_diag = 0.0;
+        int p = 0;
+        int q = 1;
+
+        for (int i = 0; i < Ndim; ++i)
+        {
+            for (int j = i + 1; j < Ndim; ++j)
+            {
+                if (std::abs(m(i, j)) > max_off_diag)
+                {
+                    max_off_diag = std::abs(m(i, j));
+                    p = i;
+                    q = j;
+                }
+            }
+        }
+
+        if (max_off_diag < eps)  break;
+
+        double phi = 0.5 * std::atan2(
+                2.0 * m(p, q), m(q, q) - m(p, p));
+        double c = std::cos(phi);
+        double s = std::sin(phi);
+
+        double m_pp = m(p, p);
+        double m_qq = m(q, q);
+        double m_pq = m(p, q);
+
+        m(p, p) = c * c * m_pp - 2.0 * s * c * m_pq + s * s * m_qq;
+        m(q, q) = s * s * m_pp + 2.0 * s * c * m_pq + c * c * m_qq;
+        m(p, q) = 0.0;
+        m(q, p) = 0.0;
+
+        for (int i = 0; i < Ndim; ++i)
+        {
+            if (i == p || i == q)  continue;
+            double m_ip = m(i, p);
+            double m_iq = m(i, q);
+            m(i, p) = c * m_ip - s * m_iq;
+            m(p, i) = m(i, p);
+            m(i, q) = s * m_ip + c * m_iq;
+            m(q, i) = m(i, q);
+        }
+
+        for (int i = 0; i < Ndim; ++i)
+        {
+            double v_ip = V(i, p);
+            double v_iq = V(i, q);
+            V(i, p) = c * v_ip - s * v_iq;
+            V(i, q) = s * v_ip + c * v_iq;
+        }
+    }
+
+    w[0] = m(0, 0);
+    w[1] = m(1, 1);
+    w[2] = m(2, 2);
+
+    for (int i = 0; i < Ndim - 1; ++i)
+    {
+        for (int j = 0; j < Ndim - 1 - i; ++j)
+        {
+            if (w[j] <= w[j + 1])  continue;
+            std::swap(w[j], w[j + 1]);
+            for (int k = 0; k < Ndim; ++k)
+            {
+                std::swap(V(k, j), V(k, j + 1));
+            }
+        }
+    }
+}
+
+
 size_t hash_value(const Vector& v)
 {
     return boost::hash_range(v.begin(), v.end());

diff --git a/src/diffpy/srreal/R3linalg.hpp b/src/diffpy/srreal/R3linalg.hpp
@@ -150,6 +150,8 @@ const Matrix& identity();
 const Matrix& zeromatrix();
 double determinant(const Matrix& A);
 const Matrix& inverse(const Matrix& A);
+void eigen_solve_3x3(const Matrix& A, Vector& eigenvalues,
+        Matrix& eigenvectors);
 
 const Vector& floor(const Vector&);
 template <class V> double norm(const V&);