From 584e8d60c9fb0e9e6edea150c97178836e4bf2e8 Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Tue, 27 Nov 2012 14:33:18 +0100
Subject: [PATCH] Calculate form factor by using an approximation of the
integral in case of large bin sizes
---
App/src/FitSuiteHelper.cpp | 2 +-
App/src/TestMesoCrystal1.cpp | 7 +-
Core/Core.pro | 3 +-
Core/FormFactors/inc/FormFactorBigCylinder.h | 61 ----------
Core/FormFactors/inc/FormFactorCrystal.h | 7 ++
.../inc/FormFactorSphereGaussianRadius.h | 4 +-
Core/FormFactors/inc/FormFactorWeighted.h | 3 +
Core/FormFactors/inc/FormFactors.h | 1 -
Core/FormFactors/inc/IFormFactor.h | 107 ++++++++++++++++--
.../FormFactors/src/FormFactorBigCylinder.cpp | 86 --------------
Core/FormFactors/src/FormFactorCylinder.cpp | 7 +-
Core/FormFactors/src/FormFactorFullSphere.cpp | 7 +-
Core/FormFactors/src/FormFactorGauss.cpp | 12 +-
Core/FormFactors/src/FormFactorLorentz.cpp | 7 +-
.../src/FormFactorParallelepiped.cpp | 4 +-
Core/FormFactors/src/FormFactorPrism3.cpp | 7 +-
Core/FormFactors/src/FormFactorPyramid.cpp | 8 +-
Core/FormFactors/src/FormFactorWeighted.cpp | 7 ++
Core/FormFactors/src/IFormFactor.cpp | 19 ++++
Core/Tools/inc/MathFunctions.h | 12 +-
20 files changed, 175 insertions(+), 196 deletions(-)
delete mode 100644 Core/FormFactors/inc/FormFactorBigCylinder.h
delete mode 100644 Core/FormFactors/src/FormFactorBigCylinder.cpp
create mode 100644 Core/FormFactors/src/IFormFactor.cpp
diff --git a/App/src/FitSuiteHelper.cpp b/App/src/FitSuiteHelper.cpp
index 2ea5dc01263..76209f477b8 100644
--- a/App/src/FitSuiteHelper.cpp
+++ b/App/src/FitSuiteHelper.cpp
@@ -117,7 +117,7 @@ void FitSuiteObserverDraw::update(IObservable *subject)
data2draw.push_back( fitObject->getChiSquaredModule()->createChi2DifferenceMap() );
// drawing
- double maximum_of_real_signal(0);
+// double maximum_of_real_signal(0);
for(size_t i_hist=0; i_hist<data2draw.size(); ++i_hist) {
c1->cd(i_hist+1);
gPad->SetLogz();
diff --git a/App/src/TestMesoCrystal1.cpp b/App/src/TestMesoCrystal1.cpp
index b8365a93141..4f6d66e1bff 100644
--- a/App/src/TestMesoCrystal1.cpp
+++ b/App/src/TestMesoCrystal1.cpp
@@ -132,10 +132,9 @@ ISample* MesoCrystalBuilder::buildSample() const
complex_t n_avg = std::sqrt(m_surface_filling_ratio*avg_n_squared_meso + 1.0 - m_surface_filling_ratio);
complex_t n_particle_adapted = std::sqrt(n_avg*n_avg + n_particle*n_particle - 1.0);
FormFactorCylinder ff_cyl(m_meso_height, m_meso_radius);
-// FormFactorBigCylinder ff_cyl(m_meso_height, m_meso_radius);
-// double bin_size = FormFactorBigCylinder::calculateBinSize(1.77*Units::angstrom, 0.051, 80);
-// std::cout << "Bin size: " << bin_size << std::endl;
-// ff_cyl.setBinSize(bin_size);
+ double bin_size_qy = IFormFactor::CalculateBinSize(1.77*Units::angstrom, 0.051, 80);
+ double bin_size_qz = IFormFactor::CalculateBinSize(1.77*Units::angstrom, 0.05, 80);
+ ff_cyl.setBinSizes(bin_size_qy, bin_size_qz);
FormFactorDecoratorDebyeWaller ff_meso(ff_cyl.clone(), m_sigma_meso_height*m_sigma_meso_height/2.0,
m_sigma_meso_radius*m_sigma_meso_radius/2.0);
diff --git a/Core/Core.pro b/Core/Core.pro
index 38b301de89f..b812ea39e97 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -41,7 +41,6 @@ SOURCES += \
Algorithms/src/OpticalFresnel.cpp \
Algorithms/src/ResolutionFunction2DSimple.cpp \
\
- FormFactors/src/FormFactorBigCylinder.cpp \
FormFactors/src/FormFactorDWBA.cpp \
FormFactors/src/FormFactorDWBAConstZ.cpp \
FormFactors/src/FormFactorDecoratorTransformation.cpp \
@@ -54,6 +53,7 @@ SOURCES += \
FormFactors/src/FormFactorLorentz.cpp \
FormFactors/src/FormFactorWeighted.cpp \
FormFactors/src/FormFactorCrystal.cpp \
+ FormFactors/src/IFormFactor.cpp \
\
Geometry/src/BasicVector3D.cpp \
Geometry/src/Normal3D.cpp \
@@ -157,7 +157,6 @@ HEADERS += \
Algorithms/inc/ResolutionFunction2DSimple.h \
Algorithms/inc/ThreadInfo.h \
\
- FormFactors/inc/FormFactorBigCylinder.h \
FormFactors/inc/FormFactorDWBA.h \
FormFactors/inc/FormFactorDWBAConstZ.h \
FormFactors/inc/FormFactorDecoratorDebyeWaller.h \
diff --git a/Core/FormFactors/inc/FormFactorBigCylinder.h b/Core/FormFactors/inc/FormFactorBigCylinder.h
deleted file mode 100644
index 0713f1d5301..00000000000
--- a/Core/FormFactors/inc/FormFactorBigCylinder.h
+++ /dev/null
@@ -1,61 +0,0 @@
-#ifndef FORMFACTORBIGCYLINDER_H_
-#define FORMFACTORBIGCYLINDER_H_
-// ********************************************************************
-// * The BornAgain project *
-// * Simulation of neutron and x-ray scattering at grazing incidence *
-// * *
-// * LICENSE AND DISCLAIMER *
-// * Lorem ipsum dolor sit amet, consectetur adipiscing elit. Mauris *
-// * eget quam orci. Quisque porta varius dui, quis posuere nibh *
-// * mollis quis. Mauris commodo rhoncus porttitor. *
-// ********************************************************************
-//! @file FormFactorBigCylinder.h
-//! @brief Definition of FormFactorBigCylinder class
-//! @author Scientific Computing Group at FRM II
-//! @date Nov 23, 2012
-
-#include "IFormFactor.h"
-#include "IStochasticParameter.h"
-
-
-class FormFactorBigCylinder : public IFormFactorBorn
-{
-public:
- FormFactorBigCylinder(double height, double radius);
-
- ~FormFactorBigCylinder();
- virtual FormFactorBigCylinder *clone() const;
-
- virtual int getNumberOfStochasticParameters() const { return 2; }
-
- virtual double getHeight() const { return m_height; }
-
- void setBinSize(double bin_size) { m_bin_size = bin_size; }
-
- static double calculateBinSize(double lambda, double phi_range, size_t n_phi);
-
-protected:
- virtual complex_t evaluate_for_q(const cvector_t &q) const;
-
-private:
- //! copy constructor and assignment operator are hidden since there is a clone method
- FormFactorBigCylinder(const FormFactorBigCylinder &);
- FormFactorBigCylinder &operator=(const FormFactorBigCylinder &);
-
- //! initialize pool parameters, i.e. register some of class members for later access via parameter pool
- virtual void init_parameters();
-
- //! approximation to radial function for integration
- double iTilde(double qR, void *params) const;
-
- //! print class
- void print(std::ostream &ostr) const;
-
-// StochasticParameter<double> *mp_height;
-// StochasticParameter<double> *mp_radius;
- double m_height;
- double m_radius;
- double m_bin_size;
-};
-
-#endif /* FORMFACTORBIGCYLINDER_H_ */
diff --git a/Core/FormFactors/inc/FormFactorCrystal.h b/Core/FormFactors/inc/FormFactorCrystal.h
index 22b490adb91..93338d9b3ea 100644
--- a/Core/FormFactors/inc/FormFactorCrystal.h
+++ b/Core/FormFactors/inc/FormFactorCrystal.h
@@ -33,6 +33,13 @@ public:
virtual void setAmbientRefractiveIndex(complex_t refractive_index);
+ //! propagate the bin sizes to the form factor to possibly enable large bin size approximations
+ virtual void setBinSizes(double delta_qy, double delta_qz) {
+ IFormFactor::setBinSizes(delta_qy, delta_qz);
+ mp_basis_form_factor->setBinSizes(delta_qy, delta_qz);
+ mp_meso_form_factor->setBinSizes(delta_qy, delta_qz);
+ }
+
protected:
virtual complex_t evaluate_for_q(const cvector_t &q) const;
private:
diff --git a/Core/FormFactors/inc/FormFactorSphereGaussianRadius.h b/Core/FormFactors/inc/FormFactorSphereGaussianRadius.h
index eb07a88b605..a8b83fcbecb 100644
--- a/Core/FormFactors/inc/FormFactorSphereGaussianRadius.h
+++ b/Core/FormFactors/inc/FormFactorSphereGaussianRadius.h
@@ -57,7 +57,9 @@ inline FormFactorSphereGaussianRadius::FormFactorSphereGaussianRadius(double mea
inline FormFactorSphereGaussianRadius* FormFactorSphereGaussianRadius::clone() const
{
- return new FormFactorSphereGaussianRadius(m_mean, m_sigma);
+ FormFactorSphereGaussianRadius *p_clone = new FormFactorSphereGaussianRadius(m_mean, m_sigma);
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
inline FormFactorSphereGaussianRadius::~FormFactorSphereGaussianRadius()
diff --git a/Core/FormFactors/inc/FormFactorWeighted.h b/Core/FormFactors/inc/FormFactorWeighted.h
index d2a1240bf84..e4d9967882f 100644
--- a/Core/FormFactors/inc/FormFactorWeighted.h
+++ b/Core/FormFactors/inc/FormFactorWeighted.h
@@ -30,6 +30,9 @@ public:
virtual complex_t evaluate(const cvector_t &k_i, const cvector_t &k_f, double alpha_i, double alpha_f) const;
virtual int getNumberOfStochasticParameters() const;
+
+ //! propagate the bin sizes to the form factor to possibly enable large bin size approximations
+ virtual void setBinSizes(double delta_qy, double delta_qz);
protected:
std::vector<IFormFactor *> m_form_factors;
std::vector<double> m_weights;
diff --git a/Core/FormFactors/inc/FormFactors.h b/Core/FormFactors/inc/FormFactors.h
index 1315bcdbf34..b58e7431798 100644
--- a/Core/FormFactors/inc/FormFactors.h
+++ b/Core/FormFactors/inc/FormFactors.h
@@ -14,7 +14,6 @@
//! @author Scientific Computing Group at FRM II
//! @date Jun 27, 2012
-#include "FormFactorBigCylinder.h"
#include "FormFactorPyramid.h"
#include "FormFactorCylinder.h"
#include "FormFactorPrism3.h"
diff --git a/Core/FormFactors/inc/IFormFactor.h b/Core/FormFactors/inc/IFormFactor.h
index 6a20d1027b3..fa28a520dc1 100644
--- a/Core/FormFactors/inc/IFormFactor.h
+++ b/Core/FormFactors/inc/IFormFactor.h
@@ -16,6 +16,8 @@
#include "Types.h"
#include "ISample.h"
+#include "MemberFunctionIntegrator.h"
+#include "MathFunctions.h"
//- -------------------------------------------------------------------
@@ -25,24 +27,29 @@
class IFormFactor : public ISample
{
public:
- IFormFactor(){}
+ IFormFactor() : m_use_large_bin_approximation(false), m_bin_qy(0.0), m_bin_qz(0.0) {}
virtual ~IFormFactor() {}
virtual IFormFactor *clone() const=0;
virtual void setAmbientRefractiveIndex(complex_t refractive_index) { (void)refractive_index; }
- /// calculate scattering amplitude for complex wavevectors
- /// @param k_i incoming wavevector
- /// @param k_f outgoing wavevector
+ //! calculate scattering amplitude for complex wavevectors
+ //! @param k_i incoming wavevector
+ //! @param k_f outgoing wavevector
virtual complex_t evaluate(const cvector_t &k_i, const cvector_t &k_f, double alpha_i, double alpha_f) const=0;
- /// return number of variable/stochastic parameters
+ //! return number of variable/stochastic parameters
virtual int getNumberOfStochasticParameters() const { return 0; }
+ //! propagate the bin sizes to the form factor to possibly enable large bin size approximations
+ virtual void setBinSizes(double delta_qy, double delta_qz);
+
virtual double getVolume() const;
virtual double getHeight() const;
+ virtual double getRadius() const;
+
virtual bool isDistributedFormFactor() const { return false; }
virtual void createDistributedFormFactors(std::vector<IFormFactor *> &form_factors,
@@ -52,8 +59,25 @@ public:
(void)nbr_samples;
}
+ static double CalculateBinSize(double lambda, double phi_range, size_t n_phi);
+
+protected:
+ bool m_use_large_bin_approximation; //!< indicates if large bin size approximation should be used
+ double m_bin_qy, m_bin_qz; //!< the sizes of the bins in q space
};
+inline void IFormFactor::setBinSizes(double delta_qy, double delta_qz)
+{
+ m_bin_qy = delta_qy;
+ m_bin_qz = delta_qz;
+ if (m_bin_qy > M_PI/2.0/getRadius() || m_bin_qz > M_PI/2.0/getHeight()) {
+ m_use_large_bin_approximation = true;
+ }
+ else {
+ m_use_large_bin_approximation = false;
+ }
+}
+
inline double IFormFactor::getVolume() const
{
cvector_t zero;
@@ -66,6 +90,19 @@ inline double IFormFactor::getHeight() const
return result;
}
+inline double IFormFactor::getRadius() const
+{
+ double result = std::sqrt(getVolume()/getHeight());
+ return result;
+}
+
+inline double IFormFactor::CalculateBinSize(double lambda, double phi_range,
+ size_t n_phi)
+{
+ double k = 2.0*M_PI/lambda;
+ return k*phi_range/(n_phi-1.0);
+}
+
class IFormFactorDecorator : public IFormFactor
{
public:
@@ -75,6 +112,12 @@ public:
virtual void setAmbientRefractiveIndex(complex_t refractive_index) { if (mp_form_factor) mp_form_factor->setAmbientRefractiveIndex(refractive_index); }
+ //! propagate the bin sizes to the form factor to possibly enable large bin size approximations
+ virtual void setBinSizes(double delta_qy, double delta_qz) {
+ IFormFactor::setBinSizes(delta_qy, delta_qz);
+ mp_form_factor->setBinSizes(delta_qy, delta_qz);
+ }
+
virtual double getHeight() const { return mp_form_factor->getHeight(); }
protected:
@@ -89,22 +132,70 @@ inline IFormFactorDecorator::~IFormFactorDecorator()
class IFormFactorBorn : public IFormFactor
{
public:
- IFormFactorBorn(){}
+ IFormFactorBorn() {}
virtual ~IFormFactorBorn() {}
virtual IFormFactorBorn *clone() const=0;
virtual complex_t evaluate(const cvector_t &k_i, const cvector_t &k_f, double alpha_i, double alpha_f) const;
+
+ //! evaluate scattering amplitude for large bin sizes
+ virtual complex_t evaluateForLargeBins(const cvector_t &q) const;
protected:
- /// evaluate scattering amplitude for complex wavevector
- /// @param q wavevector transfer \f$q\equiv k_i-k_f\f$
+ //! evaluate scattering amplitude for complex wavevector
+ //! @param q wavevector transfer \f$q\equiv k_i-k_f\f$
virtual complex_t evaluate_for_q(const cvector_t &q) const=0;
+
+ //! override volume getter to avoid infinite loop caused by big bin approximation
+ virtual double getVolume() const;
+private:
+ double bigRadialPart(double qR, void *params) const;
+ double bigZPart(double qH2, void *params) const;
};
inline complex_t IFormFactorBorn::evaluate(const cvector_t &k_i, const cvector_t &k_f, double alpha_i, double alpha_f) const
{
(void)alpha_i;
(void)alpha_f;
+ if (m_use_large_bin_approximation) {
+ return evaluateForLargeBins(k_i - k_f);
+ }
return evaluate_for_q(k_i - k_f);
}
+inline complex_t IFormFactorBorn::evaluateForLargeBins(const cvector_t& q) const
+{
+ // z parameters
+ complex_t qzH2 = q.z()*getHeight()/2.0;
+ double qzH2_d = std::abs(qzH2);
+ double effective_bin_size_h = m_bin_qz*getHeight();
+
+ // radial parameters
+ double qrR = std::abs(q.magxy())*getRadius();
+ double effective_bin_size_r = m_bin_qy*getRadius();
+ double qRmin = qrR - effective_bin_size_r/2.0;
+ double qRmax = qrR + effective_bin_size_r/2.0;
+
+ // phase from the height of the particle
+ complex_t z_phase = std::exp(complex_t(0.0, 1.0)*qzH2);
+
+ // modulus of the height of the particle
+ double z_average_intensity = (bigZPart(qzH2_d + effective_bin_size_h/2.0, (void *)0)
+ - bigZPart(qzH2_d - effective_bin_size_h/2.0, (void *)0))/effective_bin_size_h;
+ double z_modulus = std::sqrt(z_average_intensity);
+
+ // modulus of the radial part
+ MemberFunctionIntegrator<IFormFactorBorn>::mem_function p_mf = &IFormFactorBorn::bigRadialPart;
+ MemberFunctionIntegrator<IFormFactorBorn> integrator(p_mf, this);
+ double average_intensity = integrator.integrate(qRmin, qRmax, (void*)0)/effective_bin_size_r;
+ double radial_part = std::sqrt(average_intensity);
+
+ return getVolume()*radial_part*z_modulus*z_phase;
+}
+
+inline double IFormFactorBorn::getVolume() const
+{
+ cvector_t zero;
+ return std::abs(evaluate_for_q(zero));
+}
+
#endif // IFORMFACTOR_H
diff --git a/Core/FormFactors/src/FormFactorBigCylinder.cpp b/Core/FormFactors/src/FormFactorBigCylinder.cpp
deleted file mode 100644
index 9037ab922e0..00000000000
--- a/Core/FormFactors/src/FormFactorBigCylinder.cpp
+++ /dev/null
@@ -1,86 +0,0 @@
-#include "FormFactorBigCylinder.h"
-#include "StochasticDiracDelta.h"
-
-#include "MathFunctions.h"
-#include "MemberFunctionIntegrator.h"
-#include "Numeric.h"
-
-FormFactorBigCylinder::FormFactorBigCylinder(double height, double radius)
-: m_bin_size(1.0)
-{
- setName("FormFactorBigCylinder");
- m_height = height;
- m_radius = radius;
- init_parameters();
-}
-
-FormFactorBigCylinder::~FormFactorBigCylinder()
-{
-}
-
-/* ************************************************************************* */
-// initialize pool parameters, i.e. register some of class members for later access via parameter pool
-/* ************************************************************************* */
-void FormFactorBigCylinder::init_parameters()
-{
- getParameterPool()->clear();
- getParameterPool()->registerParameter("height", &m_height);
- getParameterPool()->registerParameter("radius", &m_radius);
- getParameterPool()->registerParameter("bin_size", &m_bin_size);
-}
-
-
-FormFactorBigCylinder* FormFactorBigCylinder::clone() const
-{
- FormFactorBigCylinder *p_clone = new FormFactorBigCylinder(m_height, m_radius);
- p_clone->setBinSize(m_bin_size);
- return p_clone;
-}
-
-double FormFactorBigCylinder::calculateBinSize(double lambda, double phi_range,
- size_t n_phi)
-{
- double k = 2.0*M_PI/lambda;
- return k*phi_range/(n_phi-1.0);
-}
-
-complex_t FormFactorBigCylinder::evaluate_for_q(const cvector_t &q) const
-{
- double R = m_radius;
- double H = m_height;
-
- complex_t qzH_half = q.z()*H/2.0;
- complex_t z_part = H*MathFunctions::Sinc(qzH_half)*std::exp(complex_t(0.0, 1.0)*qzH_half);
-
- double qrR = std::abs(q.magxy())*R;
- double effective_bin_size = m_bin_size*R;
- double qRmin = qrR - effective_bin_size/2.0;
- double qRmax = qrR + effective_bin_size/2.0;
-
- MemberFunctionIntegrator<FormFactorBigCylinder>::mem_function p_mf = &FormFactorBigCylinder::iTilde;
- MemberFunctionIntegrator<FormFactorBigCylinder> integrator(p_mf, this);
- double average_intensity = integrator.integrate(qRmin, qRmax, (void*)0);
-
- double J1_qrR_div_qrR = std::sqrt(average_intensity);
- double radial_part = 2.0*M_PI*R*R*J1_qrR_div_qrR;
-
- return radial_part*z_part;
-}
-
-double FormFactorBigCylinder::iTilde(double qR, void *params) const
-{
- (void)params;
- static double a = 1.0/4.0;
- static double b = std::sqrt(M_PI/3.0/std::sqrt(3.0));
-
- return a/(1.0 + std::pow(std::abs(b*qR),3.0));
-}
-
-/* ************************************************************************* */
-// print class
-/* ************************************************************************* */
-void FormFactorBigCylinder::print(std::ostream &ostr) const
-{
- ISample::print(ostr);
-// ostr << " (height:"<<m_height << " radius:"<<m_radius << ")";
-}
diff --git a/Core/FormFactors/src/FormFactorCylinder.cpp b/Core/FormFactors/src/FormFactorCylinder.cpp
index 69ddd7094bb..15eb3fb2148 100644
--- a/Core/FormFactors/src/FormFactorCylinder.cpp
+++ b/Core/FormFactors/src/FormFactorCylinder.cpp
@@ -40,14 +40,13 @@ void FormFactorCylinder::init_parameters()
FormFactorCylinder* FormFactorCylinder::clone() const
{
- return new FormFactorCylinder(m_height, m_radius);
-// return new FormFactorCylinder(mp_height->clone(), mp_radius->clone());
+ FormFactorCylinder *p_clone = new FormFactorCylinder(m_height, m_radius);
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
complex_t FormFactorCylinder::evaluate_for_q(const cvector_t &q) const
{
-// double R = mp_radius->getCurrent();
-// double H = mp_height->getCurrent();
double R = m_radius;
double H = m_height;
diff --git a/Core/FormFactors/src/FormFactorFullSphere.cpp b/Core/FormFactors/src/FormFactorFullSphere.cpp
index 150168d558d..71a9a3b5f1b 100644
--- a/Core/FormFactors/src/FormFactorFullSphere.cpp
+++ b/Core/FormFactors/src/FormFactorFullSphere.cpp
@@ -20,7 +20,6 @@ FormFactorFullSphere::FormFactorFullSphere(double radius)
FormFactorFullSphere::~FormFactorFullSphere()
{
-// delete mp_radius;
}
/* ************************************************************************* */
@@ -35,13 +34,13 @@ void FormFactorFullSphere::init_parameters()
FormFactorFullSphere* FormFactorFullSphere::clone() const
{
- return new FormFactorFullSphere(m_radius);
-// return new FormFactorFullSphere(mp_radius->clone());
+ FormFactorFullSphere *p_clone = new FormFactorFullSphere(m_radius);
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
complex_t FormFactorFullSphere::evaluate_for_q(const cvector_t &q) const
{
-// double R = mp_radius->getCurrent();
complex_t qz = q.z();
double R = m_radius;
diff --git a/Core/FormFactors/src/FormFactorGauss.cpp b/Core/FormFactors/src/FormFactorGauss.cpp
index a1db00a1818..6bc841db1b0 100644
--- a/Core/FormFactors/src/FormFactorGauss.cpp
+++ b/Core/FormFactors/src/FormFactorGauss.cpp
@@ -20,8 +20,6 @@ FormFactorGauss::FormFactorGauss(double height, double width)
m_height = height;
m_width = width;
init_parameters();
- // mp_height = new StochasticDiracDelta<double>(height);
- // mp_width = new StochasticDiracDelta<double>(width);
}
//FormFactorGauss::FormFactorGauss(StochasticParameter<double> *p_height, StochasticParameter<double> *p_width)
@@ -32,8 +30,6 @@ FormFactorGauss::FormFactorGauss(double height, double width)
FormFactorGauss::~FormFactorGauss()
{
-// delete mp_height;
-// delete mp_width;
}
@@ -50,18 +46,16 @@ void FormFactorGauss::init_parameters()
FormFactorGauss* FormFactorGauss::clone() const
{
- return new FormFactorGauss(m_height, m_width);
-// return new FormFactorGauss(mp_height->clone(), mp_radius->clone());
+ FormFactorGauss *p_clone = new FormFactorGauss(m_height, m_width);
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
complex_t FormFactorGauss::evaluate_for_q(const cvector_t &q) const
{
-// double R = mp_width->getCurrent();
-// double H = mp_height->getCurrent();
double R = m_width;
double H = m_height;
-// complex_t z_phase = std::exp(complex_t(0.0, 1.0)*q.z()*H/2.0);
complex_t z_part = H*std::exp(-q.z()*q.z()*H*H/8.0/M_PI);
complex_t radial_part = R*R*std::exp(-(q.x()*q.x()+q.y()*q.y())*R*R/8.0/M_PI);
diff --git a/Core/FormFactors/src/FormFactorLorentz.cpp b/Core/FormFactors/src/FormFactorLorentz.cpp
index 2770146df3f..fc5c8eb36ae 100644
--- a/Core/FormFactors/src/FormFactorLorentz.cpp
+++ b/Core/FormFactors/src/FormFactorLorentz.cpp
@@ -50,14 +50,13 @@ void FormFactorLorentz::init_parameters()
FormFactorLorentz* FormFactorLorentz::clone() const
{
- return new FormFactorLorentz(m_height, m_width);
-// return new FormFactorLorentz(mp_height->clone(), mp_radius->clone());
+ FormFactorLorentz *p_clone = new FormFactorLorentz(m_height, m_width);
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
complex_t FormFactorLorentz::evaluate_for_q(const cvector_t &q) const
{
-// double R = mp_radius->getCurrent();
-// double H = mp_height->getCurrent();
static const double sigma2 = 4.0*std::pow(M_PI, 2.0/3.0);
double R = m_width;
double H = m_height;
diff --git a/Core/FormFactors/src/FormFactorParallelepiped.cpp b/Core/FormFactors/src/FormFactorParallelepiped.cpp
index a16fe912a6a..3e36f063d91 100644
--- a/Core/FormFactors/src/FormFactorParallelepiped.cpp
+++ b/Core/FormFactors/src/FormFactorParallelepiped.cpp
@@ -13,7 +13,9 @@ FormFactorParallelepiped::~FormFactorParallelepiped()
FormFactorParallelepiped* FormFactorParallelepiped::clone() const
{
- return new FormFactorParallelepiped(m_height, m_radius);
+ FormFactorParallelepiped *p_clone = new FormFactorParallelepiped(m_height, m_radius);
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
complex_t FormFactorParallelepiped::evaluate_for_q(const cvector_t& q) const
diff --git a/Core/FormFactors/src/FormFactorPrism3.cpp b/Core/FormFactors/src/FormFactorPrism3.cpp
index 3974e75af22..6e90ac5a705 100644
--- a/Core/FormFactors/src/FormFactorPrism3.cpp
+++ b/Core/FormFactors/src/FormFactorPrism3.cpp
@@ -41,14 +41,13 @@ void FormFactorPrism3::init_parameters()
FormFactorPrism3* FormFactorPrism3::clone() const
{
- return new FormFactorPrism3(m_height, m_half_side );
-// return new FormFactorPrism3(mp_height->clone(), mp_half_side->clone());
+ FormFactorPrism3 *p_clone = new FormFactorPrism3(m_height, m_half_side );
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
complex_t FormFactorPrism3::evaluate_for_q(const cvector_t &q) const
{
-// double R = mp_half_side->getCurrent();
-// double H = mp_height->getCurrent();
complex_t qz = q.z();
double R = m_half_side;
double H = m_height;
diff --git a/Core/FormFactors/src/FormFactorPyramid.cpp b/Core/FormFactors/src/FormFactorPyramid.cpp
index 2359d39cd8a..1914cd6034d 100644
--- a/Core/FormFactors/src/FormFactorPyramid.cpp
+++ b/Core/FormFactors/src/FormFactorPyramid.cpp
@@ -47,15 +47,13 @@ void FormFactorPyramid::init_parameters()
FormFactorPyramid* FormFactorPyramid::clone() const
{
- return new FormFactorPyramid(m_height, m_half_side, m_alpha);
-// return new FormFactorPyramid(mp_height->clone(), mp_half_side->clone(), mp_alpha->clone());
+ FormFactorPyramid *p_clone = new FormFactorPyramid(m_height, m_half_side, m_alpha);
+ p_clone->setBinSizes(m_bin_qy, m_bin_qz);
+ return p_clone;
}
complex_t FormFactorPyramid::evaluate_for_q(const cvector_t &q) const
{
-// double H = mp_height->getCurrent();
-// double R = mp_half_side->getCurrent();
-// double tga = std::tan(mp_alpha->getCurrent());
double H = m_height;
double R = m_half_side;
diff --git a/Core/FormFactors/src/FormFactorWeighted.cpp b/Core/FormFactors/src/FormFactorWeighted.cpp
index ad403301aef..5ad7499d63f 100644
--- a/Core/FormFactors/src/FormFactorWeighted.cpp
+++ b/Core/FormFactors/src/FormFactorWeighted.cpp
@@ -53,3 +53,10 @@ int FormFactorWeighted::getNumberOfStochasticParameters() const
}
return result;
}
+
+void FormFactorWeighted::setBinSizes(double delta_qy, double delta_qz)
+{
+ for (size_t index=0; index<m_form_factors.size(); ++index) {
+ m_form_factors[index]->setBinSizes(delta_qy, delta_qz);
+ }
+}
diff --git a/Core/FormFactors/src/IFormFactor.cpp b/Core/FormFactors/src/IFormFactor.cpp
new file mode 100644
index 00000000000..3c9c5dbfc82
--- /dev/null
+++ b/Core/FormFactors/src/IFormFactor.cpp
@@ -0,0 +1,19 @@
+#include "IFormFactor.h"
+
+double IFormFactorBorn::bigRadialPart(double qR, void *params) const
+{
+ (void)params;
+ static double a = 1.0/4.0;
+ static double b = std::sqrt(M_PI/3.0/std::sqrt(3.0));
+
+ return a/(1.0 + std::pow(std::abs(b*qR),3.0));
+}
+
+double IFormFactorBorn::bigZPart(double qH2, void* params) const
+{
+ (void)params;
+ if (qH2<Numeric::double_epsilon) return qH2;
+ double qH = 2.0*qH2;
+ return (qH*MathFunctions::Si(qH) + std::cos(qH) - 1.0)/qH;
+}
+
diff --git a/Core/Tools/inc/MathFunctions.h b/Core/Tools/inc/MathFunctions.h
index 83011e34a2f..ab376bdc67a 100644
--- a/Core/Tools/inc/MathFunctions.h
+++ b/Core/Tools/inc/MathFunctions.h
@@ -23,6 +23,7 @@
#include "gsl/gsl_sf_bessel.h"
#include "gsl/gsl_sf_trig.h"
+#include "gsl/gsl_sf_expint.h"
#include "gsl/gsl_integration.h"
namespace MathFunctions
@@ -39,12 +40,16 @@ double GenerateStandardNormalRandom();
double GenerateUniformRandom();
+//! Bessel function of the first kind and order 1
double Bessel_J1(double value);
-//complex_t Bessel_J1(complex_t value);
+//! Sine integral function: \f$Si(x)\equiv\int_0^x du \sin(u)/u\f$
+double Si(double value);
+//! Sinc function: \f$Sinc(x)\equiv\sin(x)/x\f$
double Sinc(double value);
+//! Complex Sinc function: \f$Sinc(x)\equiv\sin(x)/x\f$
complex_t Sinc(complex_t value);
complex_t Laue(complex_t value, size_t N);
@@ -92,6 +97,11 @@ inline double MathFunctions::Bessel_J1(double value)
return gsl_sf_bessel_J1(value);
}
+inline double MathFunctions::Si(double value) // int_0^x du Sin(u)/u
+{
+ return gsl_sf_Si(value);
+}
+
inline double MathFunctions::Sinc(double value) // Sin(x)/x
{
return gsl_sf_sinc(value/M_PI);
--
GitLab