From 3485f0fd068ecd46c67e86f012ace2925b82af84 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (h)" <j.wuttke@fz-juelich.de>
Date: Fri, 6 May 2016 16:59:44 +0200
Subject: [PATCH] + ff triangle, for test purposes; corr and improved runff.
---
Core/FormFactors/FormFactorPolyhedron.cpp | 40 +++++++++++------
Core/FormFactors/FormFactorPolyhedron.h | 17 ++++++++
Core/FormFactors/FormFactorTriangle.cpp | 52 +++++++++++++++++++++++
Core/FormFactors/FormFactorTriangle.h | 39 +++++++++++++++++
dev-tools/math/fftest/runff.cpp | 7 ++-
5 files changed, 140 insertions(+), 15 deletions(-)
create mode 100644 Core/FormFactors/FormFactorTriangle.cpp
create mode 100644 Core/FormFactors/FormFactorTriangle.h
diff --git a/Core/FormFactors/FormFactorPolyhedron.cpp b/Core/FormFactors/FormFactorPolyhedron.cpp
index 1ce17009b4c..acacdfaa243 100644
--- a/Core/FormFactors/FormFactorPolyhedron.cpp
+++ b/Core/FormFactors/FormFactorPolyhedron.cpp
@@ -228,21 +228,17 @@ complex_t PolyhedralFace::ff_n( int n, const cvector_t q ) const
}
}
-//! Returns the contribution of this face to the polyhedral form factor.
+//! Returns the contribution ff(q) of this face to the polyhedral form factor.
complex_t PolyhedralFace::ff( const cvector_t q, const bool sym_Ci ) const
{
- //std::cout<<"Face "<<m_normal<<":\n";
- complex_t qn = q.dot(m_normal); // conj(q)*normal (BasicVector3D::dot is antilinear in 'this' argument)
- if ( std::abs(qn)<eps*q.mag() )
- return 0;
complex_t qperp;
cvector_t qpa;
decompose_q( q, qperp, qpa );
double qpa_red = m_radius_2d * qpa.mag();
complex_t qr_perp = qperp*m_rperp;
if ( qpa_red==0 ) {
- return qn * (sym_Ci ? 2.*I*sin(qr_perp) : exp(I*qr_perp)) * m_area;
+ return (sym_Ci ? 2.*I*sin(qr_perp) : exp(I*qr_perp)) * m_area;
} else if ( qpa_red < qpa_limit_series && !sym_S2 ) {
// summation of power series
#ifdef POLYHEDRAL_DIAGNOSTIC
@@ -251,10 +247,10 @@ complex_t PolyhedralFace::ff( const cvector_t q, const bool sym_Ci ) const
complex_t fac_even;
complex_t fac_odd;
if( sym_Ci ) {
- fac_even = qn * 2. * I * sin(qr_perp);
- fac_odd = qn * 2. * cos(qr_perp);
+ fac_even = 2. * I * sin(qr_perp);
+ fac_odd = 2. * cos(qr_perp);
} else {
- fac_even = qn * exp( I*qr_perp );
+ fac_even = exp( I*qr_perp );
fac_odd = fac_even;
}
complex_t sum = fac_even * m_area;
@@ -277,8 +273,8 @@ complex_t PolyhedralFace::ff( const cvector_t q, const bool sym_Ci ) const
#endif
} else {
// direct evaluation of analytic formula
- cvector_t prevec = 2.*m_normal.cross( qpa ); // complex conjugation will take place in .dot
- complex_t prefac = qn;
+ cvector_t prevec = m_normal.cross( qpa ); // complex conjugation will take place in .dot
+ complex_t prefac = 2.;
if( sym_S2 )
prefac *= sym_Ci ? -4.*sin(qr_perp) : 2.*I*exp(I*qr_perp);
complex_t sum = 0;
@@ -288,6 +284,7 @@ complex_t PolyhedralFace::ff( const cvector_t q, const bool sym_Ci ) const
complex_t Rfac = sym_S2 ? sin(e.qR(qpa)) : ( sym_Ci ? 2.*cos(qR) : exp(I*qR) );
sum += prevec.dot(e.E()) * MathFunctions::sinc(qE) * Rfac;
}
+ std::cerr<<std::setprecision(12)<<"DEBUG "<<qpa<<" "<<prevec.mag()<<" "<<sum.real()<<" "<<sum.imag()<<"\n";
//std::cout<<std::setprecision(16)<<" ret="<<prefac * sum / ( I*qpa.mag2() )<<"\n";
return prefac * sum / ( I*qpa.mag2() );
}
@@ -447,8 +444,12 @@ complex_t FormFactorPolyhedron::evaluate_centered( const cvector_t q ) const
} else {
// direct evaluation of analytic formula (coefficients may involve series)
complex_t sum = 0;
- for( const PolyhedralFace& Gk: m_faces )
- sum += Gk.ff(q, m_sym_Ci );
+ for( const PolyhedralFace& Gk: m_faces ) {
+ complex_t qn = Gk.normalProjectionConj( q ); // conj(q)*normal
+ if ( std::abs(qn)<eps*q.mag() )
+ continue;
+ sum += qn * Gk.ff(q, m_sym_Ci );
+ }
return sum / (I * q.mag2());
}
}
@@ -495,3 +496,16 @@ complex_t FormFactorPolygonalPrism::evaluate_for_q( const cvector_t q ) const
return m_height * exp(I*(m_height/2)*q.z()) * MathFunctions::sinc(m_height/2*q.z()) *
m_base->ff_2D( qxy );
}
+
+
+//***************************************************************************************************
+// FormFactorPolygonalSurface implementation
+//***************************************************************************************************
+
+complex_t FormFactorPolygonalSurface::evaluate_for_q( const cvector_t q ) const
+{
+#ifdef POLYHEDRAL_DIAGNOSTIC
+ diagnosis = { 0, 0 };
+#endif
+ return m_base->ff( q, false );
+}
diff --git a/Core/FormFactors/FormFactorPolyhedron.h b/Core/FormFactors/FormFactorPolyhedron.h
index f0f41253343..9cae7cbf56c 100644
--- a/Core/FormFactors/FormFactorPolyhedron.h
+++ b/Core/FormFactors/FormFactorPolyhedron.h
@@ -60,6 +60,8 @@ public:
kvector_t center() const { return m_center; }
double pyramidalVolume() const { return m_rperp*m_area/3; }
double radius3d() const { return m_radius_3d; }
+ //! Returns conj(q)*normal [BasicVector3D::dot is antilinear in 'this' argument]
+ complex_t normalProjectionConj( cvector_t q ) const { return q.dot(m_normal); }
complex_t ff_n( int m, const cvector_t q ) const;
complex_t ff( const cvector_t q, const bool sym_Ci ) const;
complex_t ff_2D( const cvector_t qpa ) const;
@@ -141,4 +143,19 @@ protected:
double m_height;
};
+
+//! A polygonal surface, for testing form factor computations.
+
+class FormFactorPolygonalSurface : public IFormFactorBorn {
+public:
+ FormFactorPolygonalSurface() {}
+
+ virtual complex_t evaluate_for_q(const cvector_t q ) const final;
+ double getVolume() const { return 0; }
+ virtual double getRadius() const final { return std::sqrt(m_base->area()); }
+
+protected:
+ std::unique_ptr<PolyhedralFace> m_base;
+};
+
#endif // FORMFACTORPOLYHEDRON_H
diff --git a/Core/FormFactors/FormFactorTriangle.cpp b/Core/FormFactors/FormFactorTriangle.cpp
new file mode 100644
index 00000000000..79433812741
--- /dev/null
+++ b/Core/FormFactors/FormFactorTriangle.cpp
@@ -0,0 +1,52 @@
+// ************************************************************************** //
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file FormFactors/FormFactorTriangle.cpp
+//! @brief Implements class FormFactorTriangle.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @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
+//
+// ************************************************************************** //
+
+#include "FormFactorTriangle.h"
+#include "BornAgainNamespace.h"
+#include "MathFunctions.h"
+
+//! @brief Triangle constructor
+//! @param base_edge of hexagonal base
+//! @param height of Triangle
+FormFactorTriangle::FormFactorTriangle(const double base_edge)
+ : m_base_edge( base_edge )
+{
+ setName("Triangle");
+ registerParameter(BornAgain::BaseEdge, &m_base_edge, AttLimits::n_positive());
+ onChange();
+}
+
+void FormFactorTriangle::onChange()
+{
+ double a = m_base_edge;
+ double as = a/2;
+ double ac = a/sqrt(3)/2;
+ double ah = a/sqrt(3);
+ kvector_t V[3] = {
+ { -as, -ac, 0. },
+ { as, -ac, 0. },
+ { 0., ah, 0. } };
+ m_base = std::unique_ptr<PolyhedralFace>( new PolyhedralFace( { V[0], V[1], V[2] }, false ) );
+}
+
+FormFactorTriangle* FormFactorTriangle::clone() const
+{
+ return new FormFactorTriangle(m_base_edge);
+}
+
+void FormFactorTriangle::accept(ISampleVisitor *visitor) const
+{
+ visitor->visit(this);
+}
diff --git a/Core/FormFactors/FormFactorTriangle.h b/Core/FormFactors/FormFactorTriangle.h
new file mode 100644
index 00000000000..308f6661eb4
--- /dev/null
+++ b/Core/FormFactors/FormFactorTriangle.h
@@ -0,0 +1,39 @@
+// ************************************************************************** //
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file FormFactors/FormFactorTriangle.h
+//! @brief Declares class FormFactorTriangle.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @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
+//
+// ************************************************************************** //
+
+#ifndef FORMFACTORTRIANGLE_H
+#define FORMFACTORTRIANGLE_H
+#include "FormFactorPolyhedron.h"
+
+//! @class FormFactorTriangle
+//! @ingroup formfactors
+//! @brief The formfactor of an equilateral triangle, for testing purposes.
+
+class BA_CORE_API_ FormFactorTriangle : public FormFactorPolygonalSurface
+{
+public:
+ FormFactorTriangle(const double base_edge);
+
+ virtual FormFactorTriangle *clone() const;
+ virtual void accept(ISampleVisitor *visitor) const;
+
+ double getBaseEdge() const { return m_base_edge; }
+
+private:
+ virtual void onChange() final;
+ double m_base_edge;
+};
+
+#endif // FORMFACTORTRIANGLE_H
diff --git a/dev-tools/math/fftest/runff.cpp b/dev-tools/math/fftest/runff.cpp
index e5df7d9371c..96a1418e04f 100644
--- a/dev-tools/math/fftest/runff.cpp
+++ b/dev-tools/math/fftest/runff.cpp
@@ -15,6 +15,7 @@
#include <vector>
#include "ParticleShapes.h"
+#include "FormFactorTriangle.h"
using std::cout;
using std::cerr;
@@ -61,6 +62,8 @@ IFormFactorBorn* make_particle( int ishape )
} else if( ishape==11 ) {
double alpha = 72 * Units::degree;
return new FormFactorCuboctahedron(1., 1., .8, alpha);
+ } else if( ishape==90 ) {
+ return new FormFactorTriangle(1.);
} else
throw "Shape not implemented";
}
@@ -315,10 +318,10 @@ int main (int argc, const char *argv[])
run( P, ishape, mag*uq, outfilter );
} else if( inmode==1 ) {
NEXTARG;
- mag = atof(argv[10]);
+ mag = atof( *arg );
run( P, ishape, mag*uq, outfilter );
} else if( inmode==2 ) {
- int n_mag = 2001;
+ int n_mag = 301;
double mag_i = 1e-9;
double mag_f = 1e2;
for( int i=1; i<n_mag; ++i ) {
--
GitLab