Add all relevant reciprocal peaks for angular disordered 3d lattices and...

Add all relevant reciprocal peaks for angular disordered 3d lattices and stabilize Kent distribution numerically

Core/Aggregate/IPeakShape.cpp

@@ -15,6 +15,12 @@
 #include "IPeakShape.h"
 #include "MathConstants.h"
 
+#include <limits>
+
+namespace {
+double KentDistribution(double x, double kappa);
+}
+
 IPeakShape::~IPeakShape() =default;
 
 IsotropicGaussPeakShape::IsotropicGaussPeakShape(double max_intensity, double domainsize)
@@ -86,13 +92,22 @@ double GaussKentPeakShape::evaluate(const kvector_t q, const kvector_t q_lattice
     double radial_part = m_max_intensity * std::exp(-dq2*m_radial_size*m_radial_size/2.0);
     double angular_part = 1.0;
     if (q_r*q_lat_r > 0.0) {
-        if (m_kappa>0.0) {
-            double dot_norm = q.dot(q_lattice_point)/q_r/q_lat_r;
-            double prefactor = m_kappa/(4.0*M_PI*std::sinh(m_kappa));
-            angular_part = prefactor * std::exp(m_kappa*dot_norm);
-        } else {
-            angular_part = 1.0/(4.0*M_PI);
-        }
+        double dot_norm = q.dot(q_lattice_point)/q_r/q_lat_r;
+        angular_part = KentDistribution(dot_norm, m_kappa);
     }
     return radial_part*angular_part;
 }
+
+namespace {
+double KentDistribution(double x, double kappa) {
+    static double maxkappa = std::log(1.0/std::numeric_limits<double>::epsilon())/2.0;
+    if (kappa<=0.0) {
+        return 1.0/(4.0*M_PI);
+    }
+    double prefactor = kappa / (4.0*M_PI);
+    if (kappa>maxkappa) {
+        return 2.0*prefactor*std::exp(kappa*(x-1.0));
+    }
+    return prefactor*std::exp(kappa*x)/std::sinh(kappa);
+}
+}

Core/Aggregate/InterferenceFunction3DLattice.cpp

@@ -59,11 +59,19 @@ double InterferenceFunction3DLattice::evaluate(const kvector_t q) const
     if (!mP_peak_shape)
         throw std::runtime_error("InterferenceFunction3DLattice::evaluate: "
                                  "no peak shape defined");
+    kvector_t center = q;
     double radius = 2.1 * m_rec_radius;
-    auto rec_vectors = m_lattice.reciprocalLatticeVectorsWithinRadius(q, radius);
+    double inner_radius = 0.0;
+    if (mP_peak_shape->angularDisorder()) {
+        center = kvector_t(0.0, 0.0, 0.0);
+        inner_radius = std::max(0.0, q.mag()-radius);
+        radius += q.mag();
+    }
+    auto rec_vectors = m_lattice.reciprocalLatticeVectorsWithinRadius(center, radius);
     double result = 0.0;
     for (const auto& q_rec : rec_vectors) {
-        result += mP_peak_shape->evaluate(q, q_rec)*DebyeWallerFactor(q_rec, m_dw_length);
+        if (!(q_rec.mag()<inner_radius))
+           result += mP_peak_shape->evaluate(q, q_rec)*DebyeWallerFactor(q_rec, m_dw_length);
     }
     return result;
 }