Pyramid[4] now implemented using generic polyhedron

Core/FormFactors/FormFactorPyramid.cpp

@@ -7,7 +7,7 @@
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2015
+//! @copyright Forschungszentrum Jülich GmbH 2015-16
 //! @authors   Scientific Computing Group at MLZ Garching
 //! @authors   C. Durniak, M. Ganeva, G. Pospelov, W. Van Herck, J. Wuttke
 //
@@ -15,14 +15,11 @@
 
 #include "FormFactorPyramid.h"
 #include "BornAgainNamespace.h"
-#include "MathFunctions.h"
 
-using namespace  BornAgain;
-
-FormFactorPyramid::FormFactorPyramid(
-    double length, double height, double alpha)
+FormFactorPyramid::FormFactorPyramid(double length, double height, double alpha)
+    : FormFactorPolyhedron( polyhedral_faces( length, height, alpha ), 0. )
 {
-    setName(FFPyramidType);
+    setName(BornAgain::FFPyramidType);
     m_length = length;
     m_height = height;
     m_alpha = alpha;
@@ -30,34 +27,46 @@ FormFactorPyramid::FormFactorPyramid(
     init_parameters();
 }
 
+FormFactorPyramid::~FormFactorPyramid() {}
+
+std::vector<PolyhedralFace> FormFactorPyramid::polyhedral_faces(
+    double length, double height, double alpha)
+{
+    double a = length/2;
+    double b = length/2 - height/std::tan(alpha);
+
+    kvector_t V[8] = {
+        // base:
+        { -a, -a, 0. },
+        {  a, -a, 0. },
+        {  a,  a, 0. },
+        { -a,  a, 0. },
+        // top:
+        { -b, -b, height },
+        {  b, -b, height },
+        {  b,  b, height },
+        { -b,  b, height } };
+    std::vector<PolyhedralFace> faces;
+    faces.push_back( PolyhedralFace( { V[3], V[2], V[1], V[0] }, true ) );
+    faces.push_back( PolyhedralFace( { V[0], V[1], V[5], V[4] } ) );
+    faces.push_back( PolyhedralFace( { V[1], V[2], V[6], V[5] } ) );
+    faces.push_back( PolyhedralFace( { V[2], V[3], V[7], V[6] } ) );
+    faces.push_back( PolyhedralFace( { V[3], V[0], V[4], V[7] } ) );
+    faces.push_back( PolyhedralFace( { V[4], V[5], V[6], V[7] }, true ) );
+
+    return faces;
+}
+
 bool FormFactorPyramid::check_initialization() const
 {
     bool result(true);
-    if(m_alpha<0.0 || m_alpha>Units::PID2) {
-        std::ostringstream ostr;
-        ostr << "FormFactorPyramid() -> Error in class initialization with parameters ";
-        ostr << " length:" << m_length;
-        ostr << " height:" << m_height;
-        ostr << " alpha:" << m_alpha << "\n\n";
-        ostr << "Check for '0 <= alpha <= pi/2' failed.";
-        throw Exceptions::ClassInitializationException(ostr.str());
-    }
-    if(m_length<0.0 || m_height<0.0) {
-        std::ostringstream ostr;
-        ostr << "FormFactorPyramid() -> Error in class initialization with parameters ";
-        ostr << " length:" << m_length;
-        ostr << " height:" << m_height;
-        ostr << " alpha:" << m_alpha << "\n\n";
-        ostr << "Check for '0 <= length and 0 <= height' failed.";
-        throw Exceptions::ClassInitializationException(ostr.str());
-    }
-    if(2.*m_height > m_length*std::tan(m_alpha)) {
+    if(m_height > m_length*std::tan(m_alpha)) {
         std::ostringstream ostr;
-        ostr << "FormFactorPyramid() -> Error in class initialization with parameters ";
+        ostr << "FormFactorPyramid() -> Error in class initialization with parameters";
         ostr << " length:" << m_length;
         ostr << " height:" << m_height;
-        ostr << " alpha:" << m_alpha << "\n\n";
-        ostr << "Check for '2.*height <= length*tan(alpha)' failed.";
+        ostr << " alpha[rad]:" << m_alpha << "\n\n";
+        ostr << "Check for 'height <= length*tan(alpha)' failed.";
         throw Exceptions::ClassInitializationException(ostr.str());
     }
     return result;
@@ -66,102 +75,17 @@ bool FormFactorPyramid::check_initialization() const
 void FormFactorPyramid::init_parameters()
 {
     clearParameterPool();
-    registerParameter(Length, &m_length, AttLimits::n_positive());
-    registerParameter(Height, &m_height, AttLimits::n_positive());
-    registerParameter(Alpha, &m_alpha, AttLimits::n_positive());
+    registerParameter(BornAgain::Length, &m_length, AttLimits::n_positive());
+    registerParameter(BornAgain::Height, &m_height, AttLimits::n_positive());
+    registerParameter(BornAgain::Alpha, &m_alpha, AttLimits::n_positive());
 }
 
 FormFactorPyramid* FormFactorPyramid::clone() const
 {
-    return new FormFactorPyramid(m_length, m_height, m_alpha);
+   return new FormFactorPyramid(m_length, m_height, m_alpha);
 }
 
 void FormFactorPyramid::accept(ISampleVisitor *visitor) const
 {
     visitor->visit(this);
 }
-
-double FormFactorPyramid::getRadius() const
-{
-    return m_length/2.0;
-}
-
-complex_t FormFactorPyramid::evaluate_for_q(const cvector_t& q) const
-{
-
-    double H = m_height;
-    double R = m_length/2.;
-    double tga = std::tan(m_alpha);
-    if (m_height == 0.0) return 0.0;
-    if (m_alpha == 0.0) return 0.0;
-    if (m_length == 0.0) return 0.0;
-
-    complex_t qx = q.x();
-    complex_t qy = q.y();
-    complex_t qz = q.z();
-
-    const complex_t im(0, 1);
-    complex_t full = fullPyramidPrimitive(qx/tga, qy/tga, qz, -R*tga);
-    complex_t top = fullPyramidPrimitive(qx/tga, qy/tga, qz, H-R*tga);
-    return std::exp(im*qz*R*tga) * (full-top) / (tga*tga);
-}
-
-complex_t FormFactorPyramid::fullPyramidPrimitive(complex_t a, complex_t b, complex_t c,
-                                                  double z) const
-{
-    const complex_t im(0, 1);
-    if (std::norm(a * z) > Numeric::double_epsilon && std::norm(b * z) > Numeric::double_epsilon) {
-        if (std::abs((a - b) * (a - b) - c * c) * z * z > Numeric::double_epsilon
-            && std::abs((a + b) * (a + b) - c * c) * z * z > Numeric::double_epsilon) {
-            complex_t phase = std::exp(im * c * z);
-            complex_t numerator
-                = std::sin(a * z) * (b * (a * a - b * b + c * c) * std::cos(b * z)
-                                     + im * c * (a * a + b * b - c * c) * std::sin(b * z))
-                  + a * std::cos(a * z) * ((-a * a + b * b + c * c) * std::sin(b * z)
-                                           + 2.0 * im * b * c * std::cos(b * z));
-            complex_t denominator = a * b * (a - b - c) * (a + b - c) * (a - b + c) * (a + b + c);
-            return -4.0 * phase * numerator / denominator;
-        } else {
-            return 2.0 * (g(a - b, c, z) - g(a + b, c, z)) / (a * b);
-        }
-    } else if (std::norm(a * z) <= Numeric::double_epsilon
-               && std::norm(b * z) <= Numeric::double_epsilon) {
-        if (std::norm(c * z) <= Numeric::double_epsilon) {
-            return -4.0 * std::pow(z, 3) / 3.0;
-        } else
-            return 4.0 * im
-                   * (2.0 + std::exp(im * c * z) * (c * c * z * z + 2.0 * im * c * z - 2.0))
-                   / std::pow(c, 3);
-    } else {
-        complex_t abmax;
-        if (std::norm(b * z) <= Numeric::double_epsilon) {
-            abmax = a;
-        } else {
-            abmax = b;
-        }
-        return 2.0 * (h(c - abmax, z) - h(c + abmax, z)) / abmax;
-    }
-}
-
-complex_t FormFactorPyramid::g(complex_t x, complex_t c, double z) const
-{
-    const complex_t im(0, 1);
-    if (std::abs((x*x-c*c)*z*z) > Numeric::double_epsilon) {
-        return (im*c - std::exp(im*c*z)*(x*std::sin(x*z) + im*c*std::cos(x*z))) / (x*x-c*c);
-    } else {
-        if (std::norm(c*z) > Numeric::double_epsilon) {
-            return im * (std::exp(2.0*im*c*z) + 2.0*im*c*z - 1.0) / (4.0*c);
-        } else {
-            return -z;
-        }
-    }
-}
-
-complex_t FormFactorPyramid::h(complex_t x, double z) const
-{
-    const complex_t im(0, 1);
-    if (std::norm(x*z) > Numeric::double_epsilon) {
-        return (im - (im+x*z)*std::exp(im*x*z))/(x*x);
-    }
-    return -im*z*z/2.0;
-}

Core/FormFactors/FormFactorPyramid.h

@@ -7,7 +7,7 @@
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2015
+//! @copyright Forschungszentrum Jülich GmbH 2015-16
 //! @authors   Scientific Computing Group at MLZ Garching
 //! @authors   C. Durniak, M. Ganeva, G. Pospelov, W. Van Herck, J. Wuttke
 //
@@ -16,13 +16,12 @@
 #ifndef FORMFACTORPYRAMID_H
 #define FORMFACTORPYRAMID_H
 
-#include "IFormFactorBorn.h"
+#include "FormFactorPolyhedron.h"
 
 //! @class FormFactorPyramid
 //! @ingroup formfactors
-//! @brief The form factor of pyramid.
-
-class BA_CORE_API_ FormFactorPyramid : public IFormFactorBorn
+//! @brief The formfactor of a cone6.
+class BA_CORE_API_ FormFactorPyramid : public FormFactorPolyhedron
 {
 public:
     //! @brief Pyramid constructor
@@ -30,48 +29,33 @@ public:
     //! @param height of Pyramid
     //! @param angle in radians between base and facet
     FormFactorPyramid(double length, double height, double alpha);
+    virtual ~FormFactorPyramid();
 
-    virtual ~FormFactorPyramid() {}
-    virtual FormFactorPyramid *clone() const;
+    static std::vector<PolyhedralFace> polyhedral_faces(
+        double length, double height,  double alpha);
 
-    virtual void accept(ISampleVisitor *visitor) const;
+    virtual FormFactorPyramid* clone() const;
 
-    virtual double getRadius() const;
+    virtual void accept(ISampleVisitor *visitor) const;
 
+    virtual double getRadius() const final;
     double getHeight() const;
-
     double getLength() const;
-
     double getAlpha() const;
 
-    virtual complex_t evaluate_for_q(const cvector_t& q) const;
-
 protected:
     virtual bool check_initialization() const;
     virtual void init_parameters();
 
 private:
-    complex_t fullPyramidPrimitive(complex_t a, complex_t b, complex_t c, double z) const;
-    complex_t g(complex_t x, complex_t c, double z) const;  // helper function
-    complex_t h(complex_t x, double z) const;  // helper function
     double m_length;
     double m_height;
     double m_alpha;
 };
 
-inline double FormFactorPyramid::getHeight() const
-{
-    return m_height;
-}
-
-inline double FormFactorPyramid::getLength() const
-{
-    return m_length;
-}
-
-inline double FormFactorPyramid::getAlpha() const
-{
-    return m_alpha;
-}
+inline double FormFactorPyramid::getHeight() const { return m_height; }
+inline double FormFactorPyramid::getLength() const { return m_length; }
+inline double FormFactorPyramid::getAlpha() const { return m_alpha; }
+inline double FormFactorPyramid::getRadius() const { return m_length/2.0; }
 
 #endif // FORMFACTORPYRAMID_H

Doc/UserManual/FormFactors.tex

@@ -998,7 +998,7 @@ for four different angles~$\omega$ of rotation around the $z$ axis.}
 \label{fig:FFprism3Ex}
 \end{figure}
 
-\paragraph{References}\strut\\
+\paragraph{References and History}\strut\\
 Has been validated against the \E{Prism3} form factor of \IsGISAXS\
 \cite[Eq.~2.29]{Laz08} \cite[Eq.~219]{ReLL09}.
 Note the different parameterization $L= 2 R_{\rm{\Code{IsGISAXS}}}$.
@@ -1114,23 +1114,9 @@ They must fulfill
 \end{displaymath}
 
 \paragraph{Form factor, volume, horizontal section}\strut\\
-Notation:
-\begin{displaymath}
-  a\coloneqq q_x/\tan\alpha,\quad
-  b\coloneqq q_y/\tan\alpha,\quad
-  c\coloneqq q_z.
-\end{displaymath}
-Results:
-\begin{equation*}
-  \DS F= \frac{\exp(iq_zL\tan\alpha/2)}{\tan^2\alpha} \Big\{ I_p(H-L\tan\alpha/2) - I_p(-L\tan\alpha/2) \Big\}
-\end{equation*}
-Where $I_p(z)$ is an antiderivative for the final z--integration in the Fourier transform of the non--truncated pyramid shape function:
 \begin{equation*}
-\begin{array}{@{}l@{}l@{}}
-  \DS I_p(z) = & 4\exp(icz)\Big\{ \sin(az)\big[b(a^2-b^2+c^2)\cos(bz) + ic(a^2+b^2-c^2)\sin(bz)\big] \\
-  & +a\cos(az)\big[(-a^2+b^2+c^2)\sin(bz) + 2ibc\cos(bz) \big] \Big\} \\
-  & / \Big\{ab\big[(a+b)^2-c^2\big]\big[(a-b)^2-c^2\big]\Big\}
-\end{array}
+  F \text{~: computed algebraically,
+          using the generic polyhedron form factor~\cite{ba:ffp},}
 \end{equation*}
 \begin{equation*}
   V = \dfrac{1}{6}  L^3 \tan\alpha\left[ 1
@@ -1152,12 +1138,15 @@ computed with $L=10$~nm, $H=4.2$~nm and $\alpha=60^{\circ}$,
 for four different angles~$\omega$ of rotation around the $z$ axis.}
 \end{figure}
 
-\paragraph{References}\strut\\
+\paragraph{References and History}\strut\\
 Corresponds to \E{Pyramid} form factor of \IsGISAXS\
 \cite[Eq.~2.31]{Laz08} \cite[Eq.~221]{ReLL09},
-with different parametrization $L=2R_{\rm{\Code{IsGISXAXS}}}$,
-with correction of a sign error,
-and with a more compact form of $F(\q)$.
+except for different parametrization $L=2R_{\rm{\Code{IsGISXAXS}}}$
+and a corrected sign.
+
+Reimplemented in \BornAgain-1.6 using the generic form factor
+of a polygonal prism \cite{ba:ffp},
+to achieve numerical stability near the removable singularity at $q\to0$.
 
 
 %===============================================================================