Skip to content
Snippets Groups Projects
Commit a48dd44b authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

precompute two factors

parent 7949d382
No related branches found
No related tags found
No related merge requests found
......@@ -80,29 +80,29 @@ double FormFactorRipple1::getRadius() const
//! Integrand for complex formfactor.
complex_t FormFactorRipple1::Integrand(double Z) const
{
complex_t iqZ = I*m_q.z()*Z;
complex_t aa = std::acos(2.0*Z/m_height - 1.0);
return std::exp(iqZ)*aa*MathFunctions::sinc( aa*m_q.y()*m_width/(Units::PI2) );
return std::exp(m_az*Z)*aa*MathFunctions::sinc( m_ay*aa );
}
//! Complex formfactor.
complex_t FormFactorRipple1::evaluate_for_q(const cvector_t& q) const
{
m_q = q;
complex_t factor = m_length*MathFunctions::sinc(m_q.x()*m_length*0.5)*m_width/Units::PI;
complex_t factor = m_length*MathFunctions::sinc(q.x()*m_length*0.5)*m_width/Units::PI;
// analytical expressions for some particular cases
if ( m_q.z()==0. ) {
if( m_q.y()==0. )
if ( q.z()==0. ) {
if( q.y()==0. )
return factor*Units::PID2*m_height;
complex_t aaa = m_q.y()*m_width/(Units::PI2);
complex_t aaa = q.y()*m_width/(Units::PI2);
complex_t aaa2 = aaa*aaa;
if ( aaa2==1. )
return factor*Units::PID4*m_height;
return factor*Units::PID2*m_height*MathFunctions::sinc(m_q.y()*m_width*0.5)/(1.0-aaa2);
return factor*Units::PID2*m_height*MathFunctions::sinc(q.y()*m_width*0.5)/(1.0-aaa2);
}
// numerical integration otherwise
m_ay = q.y() * m_width / Units::PI2;
m_az = I*q.z();
complex_t integral = mP_integrator->integrate(0, m_height);
return factor*integral;
}
......@@ -35,7 +35,6 @@ public:
//! @param width of cosine cross section
//! @param height of cosine cross section
FormFactorRipple1(double length, double width, double height);
virtual ~FormFactorRipple1();
virtual FormFactorRipple1 *clone() const;
......@@ -43,11 +42,8 @@ public:
virtual void accept(ISampleVisitor *visitor) const;
virtual double getRadius() const;
double getHeight() const;
double getWidth() const;
double getLength() const;
virtual complex_t evaluate_for_q(const cvector_t& q) const;
......@@ -62,7 +58,9 @@ private:
double m_width;
double m_height;
double m_length;
mutable cvector_t m_q;
mutable complex_t m_ay;
mutable complex_t m_az;
#ifndef GCCXML_SKIP_THIS
std::unique_ptr<IntegratorComplex<FormFactorRipple1>> mP_integrator;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment