1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
|
// ************************************************************************************************
//
// libformfactor: efficient and accurate computation of scattering form factors
//
//! @file ff/Edge.cpp
//! @brief Implements class Edge
//!
//! @homepage https://jugit.fz-juelich.de/mlz/libformfactor
//! @license GNU General Public License v3 or higher (see LICENSE)
//! @copyright Forschungszentrum Jülich GmbH 2022
//! @author Joachim Wuttke, Scientific Computing Group at MLZ (see CITATION)
//
// ************************************************************************************************
#include "ff/Edge.h"
#include "ff/Factorial.h"
#include <algorithm>
#include <iomanip>
#include <stdexcept>
namespace {
constexpr auto ReciprocalFactorialArray = ff_aux::generateReciprocalFactorialArray<171>();
} // namespace
ff::Edge::Edge(R3 Vlow, R3 Vhig) : m_E((Vhig - Vlow) / 2), m_R((Vhig + Vlow) / 2) {}
//! Returns sum_l=0^M/2 u^2l v^(M-2l) / (2l+1)!(M-2l)! - vperp^M/M!
complex_t ff::Edge::contrib(int M, C3 qpa, complex_t qrperp) const
{
complex_t u = qE(qpa);
complex_t v2 = m_R.dot(qpa);
complex_t v1 = qrperp;
complex_t v = v2 + v1;
// std::cout << std::scientific << std::showpos << std::setprecision(16) << "contrib: u=" << u
// << " v1=" << v1 << " v2=" << v2 << "\n";
if (v == 0.) { // only 2l=M contributes
if (M & 1) // M is odd
return 0.;
return ReciprocalFactorialArray[M] * (pow(u, M) / (M + 1.) - pow(v1, M));
}
complex_t result = 0;
// the l=0 term, minus (qperp.R)^M, which cancels under the sum over E*contrib()
if (v1 == 0.)
result = ReciprocalFactorialArray[M] * pow(v2, M);
else if (v2 == 0.) {
; // leave result=0
} else {
// binomial expansion
for (int mm = 1; mm <= M; ++mm) {
complex_t term = ReciprocalFactorialArray[mm] * ReciprocalFactorialArray[M - mm]
* pow(v2, mm) * pow(v1, M - mm);
result += term;
// std::cout << "contrib mm=" << mm << " t=" << term << " s=" << result << "\n";
}
}
if (u == 0.)
return result;
for (int l = 1; l <= M / 2; ++l) {
complex_t term = ReciprocalFactorialArray[M - 2 * l] * ReciprocalFactorialArray[2 * l + 1]
* pow(u, 2 * l) * pow(v, M - 2 * l);
result += term;
// std::cout << "contrib l=" << l << " t=" << term << " s=" << result << "\n";
}
return result;
}
|