Adapted implementation of IF2DParaCrystal

Core/Samples/inc/InterferenceFunction2DParaCrystal.h

@@ -34,16 +34,19 @@ public:
     //! @param alpha_lattice Angle between lattice basis vectors.
     //! @param xi Angle between first basis vector and the x-axis of incoming beam.
     //! @param m_corr_length correlation length of paracrystal
-    InterferenceFunction2DParaCrystal(double length_1, double length_2, double alpha_lattice, double xi=0.0, double corr_length=0.0);
+    InterferenceFunction2DParaCrystal(double length_1, double length_2,
+            double alpha_lattice, double xi=0.0, double corr_length=0.0);
     virtual ~InterferenceFunction2DParaCrystal();
 
     virtual InterferenceFunction2DParaCrystal *clone() const;
 
     virtual void accept(ISampleVisitor *visitor) const { visitor->visit(this); }
 
-    static InterferenceFunction2DParaCrystal *createSquare(double peak_distance, double corr_length=0.0,
+    static InterferenceFunction2DParaCrystal *createSquare(
+            double peak_distance, double corr_length=0.0,
             double domain_size_1=0.0, double domain_size_2=0.0);
-    static InterferenceFunction2DParaCrystal *createHexagonal(double peak_distance, double corr_length=0.0,
+    static InterferenceFunction2DParaCrystal *createHexagonal(
+            double peak_distance, double corr_length=0.0,
             double domain_size_1=0.0, double domain_size_2=0.0);
 
     //! @brief Sets sizes of coherence domain
@@ -53,20 +56,26 @@ public:
         m_domain_sizes[1] = size_2;
     }
 
-    void setProbabilityDistributions(const IFTDistribution2D& pdf_1, const IFTDistribution2D& pdf_2);
+    void setProbabilityDistributions(const IFTDistribution2D& pdf_1,
+            const IFTDistribution2D& pdf_2);
 
-    void setIntegrationOverXi(bool integrate_xi) { m_integrate_xi = integrate_xi; }
+    void setIntegrationOverXi(bool integrate_xi)
+    { m_integrate_xi = integrate_xi; }
 
     virtual double evaluate(const cvector_t& q) const;
 
-    //! Adds parameters from local pool to external pool and call recursion over direct children
-    virtual std::string addParametersToExternalPool(std::string path, ParameterPool *external_pool, int copy_number=-1) const;
+    //! Adds parameters from local pool to external pool and call recursion
+    //! over direct children
+    virtual std::string addParametersToExternalPool(std::string path,
+            ParameterPool *external_pool, int copy_number=-1) const;
+
 
 protected:
     //! Registers some class members for later access via parameter pool
     virtual void init_parameters();
 
-    void transformToPrincipalAxes(double qx, double qy, double gamma, double delta, double& q_pa_1, double& q_pa_2) const;
+    void transformToPrincipalAxes(double qx, double qy, double gamma,
+            double delta, double& q_pa_1, double& q_pa_2) const;
     double m_lattice_lengths[2]; //!< the size of unit cell
     double m_alpha_lattice; //!< Angle between lattice basis vectors
     double m_xi; //!< Orientation of the lattice wrt beam axis x
@@ -83,6 +92,9 @@ private:
     //! Returns interference function for fixed xi in 1d
     double interference1D(double qx, double qy, double xi, size_t index) const;
 
+    //! Calculates the geometric series of z to order N
+    complex_t geometricSum(complex_t z, int exponent) const;
+
     complex_t FTPDF(double qx, double qy, double xi, size_t index) const;
 
     mutable double m_qx;

Core/Samples/src/InterferenceFunction2DParaCrystal.cpp

@@ -179,15 +179,20 @@ double InterferenceFunction2DParaCrystal::interference1D(double qx, double qy,
     if (n<1) {
         result = ((1.0 + fp)/(1.0 - fp)).real();
     } else {
-        if (std::abs(1.0-fp) < Numeric::double_epsilon) {
+        if (std::abs(1.0-fp) < Numeric::double_epsilon ) {
             result = nd;
-        } else {
+        }
+        else if (std::abs(1.0-fp)*nd < 2e-4) {
+            double intermediate = geometricSum(fp, n).real()/nd;
+            result = 1.0 + 2.0*intermediate;
+        }
+        else {
             complex_t tmp;
             double double_min = std::numeric_limits<double>::min();
             if (std::log(std::abs(fp)+double_min)*nd < std::log(double_min)) {
                 tmp = complex_t(0.0, 0.0);
             } else {
-                tmp = std::pow(fp,n-1);
+            tmp = std::pow(fp,n-1);
             }
             double intermediate = ((1.0-1.0/nd)*fp/(1.0-fp)
                     - fp*fp*(1.0-tmp)/nd/(1.0-fp)/(1.0-fp)).real();
@@ -197,6 +202,24 @@ double InterferenceFunction2DParaCrystal::interference1D(double qx, double qy,
     return result;
 }
 
+complex_t InterferenceFunction2DParaCrystal::geometricSum(complex_t z,
+        int exponent) const
+{
+    if (exponent<1) {
+        throw LogicErrorException(
+                "InterferenceFunction2DParaCrystal::geometricSeries:"
+                " exponent should be > 0");
+    }
+    complex_t result(0.0, 0.0);
+    double nd = (double)exponent;
+    --exponent;
+    while (exponent>0) {
+        result += std::pow(z, exponent)*(nd-exponent);
+        --exponent;
+    }
+    return result;
+}
+
 complex_t InterferenceFunction2DParaCrystal::FTPDF(double qx, double qy,
         double xi, size_t index) const
 {