From afcc26ce9348545dff5cc334361d72e62bb84aa5 Mon Sep 17 00:00:00 2001
From: Celine Durniak <c.durniak@fz-juelich.de>
Date: Tue, 26 Nov 2013 12:48:14 +0100
Subject: [PATCH] New form factor of cone6

---
 Core/FormFactors/inc/FormFactorCone6.h   |  69 +++++++++++++++
 Core/FormFactors/src/FormFactorCone6.cpp | 105 +++++++++++++++++++++++
 2 files changed, 174 insertions(+)
 create mode 100644 Core/FormFactors/inc/FormFactorCone6.h
 create mode 100644 Core/FormFactors/src/FormFactorCone6.cpp

diff --git a/Core/FormFactors/inc/FormFactorCone6.h b/Core/FormFactors/inc/FormFactorCone6.h
new file mode 100644
index 00000000000..02230939d78
--- /dev/null
+++ b/Core/FormFactors/inc/FormFactorCone6.h
@@ -0,0 +1,69 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      FormFactors/inc/FormFactorCone6.h
+//! @brief     Defines class FormFactorCone6
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#ifndef FORMFACTORCONE6_H
+#define FORMFACTORCONE6_H
+
+#include "IFormFactorBorn.h"
+#include "IStochasticParameter.h"
+#include "MemberComplexFunctionIntegrator.h"
+
+//! Form factor of a cone6.
+
+class BA_CORE_API_ FormFactorCone6 : public IFormFactorBorn
+{
+public:
+    //! @brief cone6 constructor
+    //! @param radius of hexagonal base (different from R in IsGisaxs)
+    //! @param height of cone6
+    //! @param angle in radians between base and facet
+    FormFactorCone6(double radius, double height,  double alpha);
+    ~FormFactorCone6() {delete m_integrator;}
+
+    virtual FormFactorCone6* clone() const;
+
+    virtual void accept(ISampleVisitor *visitor) const { visitor->visit(this); }
+
+    virtual int getNumberOfStochasticParameters() const { return 3; }
+
+    virtual double getHeight() const { return m_height; }
+    virtual void setHeight(double height) { m_height = height; }
+
+    virtual double getRadius() const { return m_radius; }
+    virtual void setRadius(double radius) { m_radius = radius; }
+
+    virtual double getAlpha() const { return m_alpha; }
+    virtual void setAlpha(double alpha) { m_alpha = alpha; }
+
+protected:
+    virtual complex_t evaluate_for_q (const cvector_t& q) const;
+    virtual void init_parameters();
+
+private:
+
+    complex_t Integrand(double Z, void* params) const;
+
+    double m_radius;
+    double m_height;
+    double m_alpha;
+    double m_root3; // Cached value of square root of 3
+    mutable cvector_t m_q;
+
+    MemberComplexFunctionIntegrator<FormFactorCone6> *m_integrator;
+};
+
+#endif // FORMFACTORCONE6_H
+
+
diff --git a/Core/FormFactors/src/FormFactorCone6.cpp b/Core/FormFactors/src/FormFactorCone6.cpp
new file mode 100644
index 00000000000..57b0ed7665a
--- /dev/null
+++ b/Core/FormFactors/src/FormFactorCone6.cpp
@@ -0,0 +1,105 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      FormFactors/src/FormFactorCone6.cpp
+//! @brief     Implements class FormFactorCone6.
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#include "FormFactorCone6.h"
+#include "StochasticDiracDelta.h"
+#include "Numeric.h"
+#include "MathFunctions.h"
+#include <cmath>
+#include "MemberFunctionIntegrator.h"
+#include "MemberComplexFunctionIntegrator.h"
+
+
+FormFactorCone6::FormFactorCone6(double radius, double height, double alpha)
+{
+    setName("FormFactorCone6");
+    m_radius = radius;
+    m_height = height;
+    m_alpha = alpha;
+    m_root3 = std::sqrt(3.0);
+    init_parameters();
+
+    MemberComplexFunctionIntegrator<FormFactorCone6>::mem_function p_mf =
+       & FormFactorCone6::Integrand;
+    m_integrator =
+        new MemberComplexFunctionIntegrator<FormFactorCone6>(p_mf, this);
+}
+
+void FormFactorCone6::init_parameters()
+{
+    clearParameterPool();
+    registerParameter("radius", &m_radius);
+    registerParameter("height", &m_height);
+    registerParameter("alpha", &m_alpha);
+}
+
+FormFactorCone6* FormFactorCone6::clone() const
+{
+   FormFactorCone6* result = new FormFactorCone6(m_radius, m_height, m_alpha);
+   result->setName(getName());
+   return result;
+}
+
+
+//! Integrand for complex formfactor.
+complex_t FormFactorCone6::Integrand(double Z, void* params) const
+{
+    (void)params;  // to avoid unused-variable warning
+    double Rz = m_radius -2.*Z/std::tan(m_alpha)/m_root3;
+    complex_t qx =m_q.x();
+    complex_t qy =m_q.y();
+    complex_t qz =m_q.z();
+
+    complex_t qxR_half = qx*Rz/2.0;
+    complex_t qyr3R_half = m_root3*qy*Rz/2.;
+
+    if (std::abs(3.0*qy*qy-qx*qx)==0.0) {
+
+        return Rz*Rz*m_root3/2.0*MathFunctions::Sinc(qyr3R_half)*(
+                    MathFunctions::Sinc(qyr3R_half)
+                    + 2.0*std::cos(qyr3R_half)
+                    );
+    } else {
+        return  (3./4.*qy*Rz*qy*Rz*
+              MathFunctions::Sinc(qxR_half) *
+              MathFunctions::Sinc(qyr3R_half) +
+             std::cos(2.0*qxR_half) - std::cos(qyr3R_half) * std::cos(qxR_half)
+             )*std::exp(complex_t(0.0, 1.0)*qz*Z);}
+}
+
+
+//! Complex formfactor.
+
+complex_t FormFactorCone6::evaluate_for_q(const cvector_t& q) const
+{   m_q = q;
+
+    if ( std::abs(q.mag()) < Numeric::double_epsilon) {
+
+        double R = m_radius;
+        double H = m_height;
+        double tga = std::tan(m_alpha);
+        double HdivRtga = 2./m_root3*H/tga/R;
+        return  3.0/4.0*tga*R*R*R*
+                (1.0 - (1.0 - HdivRtga)*(1.0 - HdivRtga)*(1.0 - HdivRtga));
+
+    }
+    else {
+    complex_t integral = m_integrator->integrate(0., m_height);
+    return 4.0*m_root3/(3.0*m_q.y()*m_q.y()-m_q.x()*m_q.x())*integral;
+    }
+}
+
+
+
-- 
GitLab