ENH: add ThreeDPDFCalculator support

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)
  {

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();

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();
 

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());

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&);