From 25e36fd3c9a08762ed58a1820c478b06d5cb3196 Mon Sep 17 00:00:00 2001
From: Celine Durniak <c.durniak@fz-juelich.de>
Date: Tue, 26 Nov 2013 12:49:40 +0100
Subject: [PATCH] New form factor of Tetrahedron
---
Core/FormFactors/inc/FormFactorTetrahedron.h | 18 ++-
.../FormFactors/src/FormFactorTetrahedron.cpp | 148 +++++++++++++-----
2 files changed, 122 insertions(+), 44 deletions(-)
diff --git a/Core/FormFactors/inc/FormFactorTetrahedron.h b/Core/FormFactors/inc/FormFactorTetrahedron.h
index 768007db82c..93acaffd1b3 100644
--- a/Core/FormFactors/inc/FormFactorTetrahedron.h
+++ b/Core/FormFactors/inc/FormFactorTetrahedron.h
@@ -19,19 +19,21 @@
#include "IFormFactorBorn.h"
#include "IStochasticParameter.h"
+//#include "MemberComplexFunctionIntegrator.h"
+
//! Form factor of tetrahedron.
class BA_CORE_API_ FormFactorTetrahedron : public IFormFactorBorn
{
public:
- //! @brief tetrahedron constructor
- //! @param height of tetrahedron
- //! @param half_side half of tetrahedron's base
+ //! @brief Tetrahedron constructor
+ //! @param height of Tetrahedron
+ //! @param half_side: half of Tetrahedron's base
//! @param angle in radians between base and facet
- FormFactorTetrahedron(double height, double half_side, double alpha);
-
+ FormFactorTetrahedron(double half_side, double height, double alpha);
~FormFactorTetrahedron() {}
+
virtual FormFactorTetrahedron *clone() const;
virtual void accept(ISampleVisitor *visitor) const { visitor->visit(this); }
@@ -53,10 +55,16 @@ protected:
virtual void init_parameters();
private:
+
+ //complex_t Integrand(double Z, void* params) const;
+
double m_height;
double m_half_side;
double m_alpha;
double m_root3; // Cached value of square root of 3
+ // mutable cvector_t m_q;
+
+ // MemberComplexFunctionIntegrator<FormFactorTetrahedron> *m_integrator;
};
#endif // FORMFACTORTETRAHEDRON_H
diff --git a/Core/FormFactors/src/FormFactorTetrahedron.cpp b/Core/FormFactors/src/FormFactorTetrahedron.cpp
index d94b557de3c..a7e245e60e8 100644
--- a/Core/FormFactors/src/FormFactorTetrahedron.cpp
+++ b/Core/FormFactors/src/FormFactorTetrahedron.cpp
@@ -18,7 +18,7 @@
#include "MathFunctions.h"
FormFactorTetrahedron::FormFactorTetrahedron(
- double height, double half_side, double alpha)
+ double half_side, double height, double alpha)
{
setName("FormFactorTetrahedron");
m_height = height;
@@ -26,6 +26,11 @@ FormFactorTetrahedron::FormFactorTetrahedron(
m_alpha = alpha;
m_root3 = std::sqrt(3.0);
init_parameters();
+
+ /* MemberComplexFunctionIntegrator<FormFactorTetrahedron>::mem_function p_mf =
+ & FormFactorTetrahedron::Integrand;
+ m_integrator =
+ new MemberComplexFunctionIntegrator<FormFactorTetrahedron>(p_mf, this);*/
}
void FormFactorTetrahedron::init_parameters()
@@ -39,64 +44,129 @@ void FormFactorTetrahedron::init_parameters()
FormFactorTetrahedron* FormFactorTetrahedron::clone() const
{
FormFactorTetrahedron *result =
- new FormFactorTetrahedron(m_height, m_half_side, m_alpha);
+ new FormFactorTetrahedron(m_half_side, m_height, m_alpha);
result->setName(getName());
return result;
}
+
+//! Integrand for complex formfactor.
+/*complex_t FormFactorTetrahedron::Integrand(double Z, void* params) const
+{
+ (void)params; // to avoid unused-variable warning
+ double Rz = m_half_side -Z*m_root3/std::tan(m_alpha);
+ complex_t r3qyR = m_root3*m_q.y()*Rz;
+ complex_t qxR = m_q.x()*Rz;
+
+ return (std::exp(complex_t(0.0, 1.0)*r3qyR) -
+ std::cos(qxR)-complex_t(0.0, 1.0)*r3qyR*
+ MathFunctions::Sinc(qxR))
+ *std::exp(complex_t(0.0, 1.0)*(m_q.z()*Z-m_q.y()*Rz/m_root3));
+
+}*/
+
+//! Complex formfactor.
+
+/*complex_t FormFactorTetrahedron::evaluate_for_q(const cvector_t& q) const
+{ m_q = q;
+ double H = m_height;
+ double R = m_half_side;
+ double tga = std::tan(m_alpha);
+
+ if ( std::abs(m_q.mag()) < Numeric::double_epsilon) {
+ double sqrt3HdivRtga = m_root3*H/R/tga;
+ return tga/3.*R*R*R*(1. -
+ (1.-sqrt3HdivRtga)*(1.-sqrt3HdivRtga)*(1.-sqrt3HdivRtga));
+ }
+ else {
+ complex_t integral = m_integrator->integrate(0., H);
+ return 2.0*m_root3/(q.x()*q.x()-3.0*q.y()*q.y())*integral;
+ }
+}*/
+
+
complex_t FormFactorTetrahedron::evaluate_for_q(const cvector_t& q) const
{
double H = m_height;
double R = m_half_side;
double tga = std::tan(m_alpha);
+ double L = 2.*tga*R/m_root3-H;
complex_t qx = q.x();
complex_t qy = q.y();
complex_t qz = q.z();
+ complex_t q1, q2, q3;
+ q1=(1./2.)*((m_root3*qx - qy)/tga - qz);
+ q2=(1./2.)*((m_root3*qx + qy)/tga + qz);
+ q3 = (qy/tga - qz/2.);
+
complex_t F;
const complex_t im(0,1);
- if (std::abs(qx)==0.0 && std::abs(qy)==0.0) {
- complex_t qzH_half = qz*H/2.0;
- //WRONG
- F = m_root3*R*R*H*std::exp(im*qzH_half)*MathFunctions::Sinc(qzH_half);
+ if (std::abs(qx) <= Numeric::double_epsilon) {
+
+ if (std::abs(qy) <= Numeric::double_epsilon) {
+
+ if (std::abs(qz) <= Numeric::double_epsilon) {
+ // qx=qy=qz=0
+ double sqrt3HdivRtga = m_root3*H/R/tga;
+ F = tga/3.*R*R*R*(1. -
+ (1.-sqrt3HdivRtga)*(1.-sqrt3HdivRtga)*(1.-sqrt3HdivRtga));
+ }
+ else {
+ //qx=qy=0 qz!=0
+ complex_t qzH_half = qz*H/2.0;
+ F = m_root3*H*std::exp(im*qzH_half)*(
+ MathFunctions::Sinc(qzH_half)*
+ (R*R-im*2.*m_root3*R/(qz*tga)-6./(qz*qz*tga*tga))
+ + m_root3*std::exp(im*qzH_half)/(qz*tga)*(
+ 2.*im*R - im*m_root3*H/tga + 2.*m_root3/(qz*tga)));
+ }
+ }
+ else {
+ // qx=0 qy!=0
+ F = 2.*H/(m_root3*qy*qy)*
+ std::exp(im*R*qz*tga/m_root3 + im*q3*L)*(
+ std::exp(-im*q2*L)*MathFunctions::Sinc(q2*H)
+ - std::exp(im*q3*L)*MathFunctions::Sinc(q3*H)) ;
+ }
}
else {
- if (std::abs(qx*qx-3.0*qy*qy)==0.0) {
- complex_t qa = 2.*qy/tga + qz/2.;
- complex_t qb = - qy/tga + qz/2.;
- F = H/2.*m_root3*std::exp(im*2.*qy*R/m_root3)*(
- - std::exp(im*qa*H-im*2.*m_root3*qy*R)*MathFunctions::Sinc(qa*H)
- + std::exp(im*qb*H)*MathFunctions::Sinc(qb*H)*
- ( 1.- 3.*qy/(tga*qb) - im*2.*m_root3*qy*R )
- + 3.*qy/(qb*tga)*std::exp(im*2.*qb*H)
- )/( qx*qx ) ;
- } else {
-
- complex_t q1, q2, q3;
- double L;
- L = 2.*tga*R/m_root3-H;
- q3 = (qy/tga - qz/2.);
-
- if (std::abs(qx)==0.0) {
-
- F = - 2.*H*MathFunctions::Sinc(q3*H)*
- std::exp(im*(q3*L+qz*R*tga/m_root3))/((m_root3*qy*qy));
-
- } else {
-
- q1=(1./2.)*((m_root3*qx-qy)/tga - qz);
- q2=(1./2.)*((m_root3*qx+qy)/tga + qz);
-
- F = -(1. + m_root3*qy/qx)*MathFunctions::Sinc(q1*H)*std::exp(im*q1*L)
- -(1. - m_root3*qy/qx)*MathFunctions::Sinc(q2*H)*std::exp(-im*q2*L)
- + 2.*MathFunctions::Sinc(q3*H)*std::exp(im*q3*L);
-
- F = F*H*std::exp(im*qz*R*tga/m_root3)*m_root3/((qx*qx-3.*qy*qy));}
+ // qx!=0
+ if (std::abs(qx*qx-3.0*qy*qy)==0.0) {
+ // qx**2= 3qy**2
+ complex_t qa = 2.0*qy/tga + qz/2.0;
+ F = H*m_root3*std::exp(im*2.0*qy*R/m_root3)*(
+ MathFunctions::Sinc(q3*H)
+ *(1.0 +m_root3*qy*(-im*2.0*R + m_root3/(q3*tga)))
+ - 3.0*qy*std::exp(-im*2.0*q3*H)/(q3*tga)
+ - std::exp(im*qa*H-im*2.0*m_root3*qy*R)*MathFunctions::Sinc(qa*H)
+ )/(2.0*qx*qx);
}
- }
-
+ else {
+ //Formula Isgisaxs
+ /* complex_t qc1 = 2./m_root3*qx;
+ complex_t qc2 = qx/m_root3 + qy;
+
+ complex_t qq1 = 2.*qc1 - qc2 - tga*qz;
+ complex_t qq2 = qc1 + qc2 + tga*qz;
+ complex_t qq3 =-qc1 +2.*qc2 - tga*qz;
+ complex_t Rt = R/m_root3;
+ complex_t Ht = H/(2.*tga);
+
+ F = -qc2*MathFunctions::Sinc(qq1*Ht)*std::exp(im*(Rt-Ht)*qq1) +
+ (qc2-qc1)*MathFunctions::Sinc(qq2*Ht)*std::exp(-im*(Rt-Ht)*qq2) +
+ qc1*MathFunctions::Sinc(qq3*Ht)*std::exp(im*(Rt-Ht)*qq3);
+ F = F*2.*H/m_root3/qc1/qc2/(qc1-qc2)*std::exp(im*Rt*tga*qz);*/
+
+ F =-(1.+m_root3*qy/qx)*MathFunctions::Sinc(q1*H)*std::exp(im*q1*L)
+ -(1.-m_root3*qy/qx)*MathFunctions::Sinc(q2*H)*std::exp(-im*q2*L)
+ +2.*MathFunctions::Sinc(q3*H)*std::exp(im*q3*L);
+
+ F = H*m_root3*std::exp(im*qz*R*tga/m_root3)/(qx*qx-3.*qy*qy)*F;
+}
+}
return F;
}
--
GitLab