Skip to content
Snippets Groups Projects
Commit 94ec0ac8 authored by Mikhail Svechnikov's avatar Mikhail Svechnikov
Browse files

simplify summation

parent aec06e31
No related branches found
No related tags found
1 merge request!2358Optimize exponents summation in "exact" mesocrystal FF algorithm (#892)
Pipeline #129737 passed
...@@ -86,25 +86,16 @@ complex_t ReMesocrystal::structureFactor(const WavevectorInfo& wavevectors) cons ...@@ -86,25 +86,16 @@ complex_t ReMesocrystal::structureFactor(const WavevectorInfo& wavevectors) cons
complex_t exp_c = exp_I(m_lattice.basisVectorC().dot(q)); complex_t exp_c = exp_I(m_lattice.basisVectorC().dot(q));
std::vector<complex_t> exponents_b(range.y() + 1); std::vector<complex_t> exponents_b(range.y() + 1);
std::vector<complex_t> exponents_c(range.z() + 1);
for (int j = 0; j <= range.y(); j++) for (int j = 0; j <= range.y(); j++)
exponents_b[j] = pow(exp_b, j + min.y()); exponents_b[j] = pow(exp_b, j + min.y());
for (int k = 0; k <= range.z(); k++)
exponents_c[k] = pow(exp_c, k + min.z());
complex_t denom_factor = Numeric::almostEqual(exp_c, 1., 1) ? 1. : 1. / (exp_c - 1.);
std::vector<complex_t> exp_c_sum_factor(range.z() + 1);
for (int k = 0; k <= range.z(); k++)
exp_c_sum_factor[k] = (pow(exp_c, k + 1) - 1.) * denom_factor;
// main summation // main summation
complex_t field_factor(0.0, 0.0); complex_t field_factor(0.0, 0.0);
if (Numeric::almostEqual(exp_c, 1., 1)) { if (Numeric::almostEqual(exp_c, 1., 1)) {
// Special case if exp_c == 1 // if exp_c == 1, the sum is even simpler
for (int i = 0; i <= range.x(); i++) { for (int i = 0; i <= range.x(); i++) {
const complex_t e_a = pow(exp_a, i + min.x()); const complex_t e_a = pow(exp_a, i + min.x());
...@@ -117,6 +108,14 @@ complex_t ReMesocrystal::structureFactor(const WavevectorInfo& wavevectors) cons ...@@ -117,6 +108,14 @@ complex_t ReMesocrystal::structureFactor(const WavevectorInfo& wavevectors) cons
} }
} }
} else { } else {
std::vector<complex_t> exponents_c(range.z() + 1);
for (int k = 0; k <= range.z(); k++)
exponents_c[k] = pow(exp_c, k + min.z());
std::vector<complex_t> exp_c_sum_factor(range.z() + 1);
for (int k = 0; k <= range.z(); k++)
exp_c_sum_factor[k] = (pow(exp_c, k + 1) - 1.) / (exp_c - 1.);
for (int i = 0; i <= range.x(); i++) { for (int i = 0; i <= range.x(); i++) {
const complex_t e_a = pow(exp_a, i + min.x()); const complex_t e_a = pow(exp_a, i + min.x());
......
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