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
67
68
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Resample/Interparticle/DecouplingApproximationStrategy.cpp
//! @brief Implements class DecouplingApproximationStrategy.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#include "Resample/Interparticle/DecouplingApproximationStrategy.h"
#include "Base/Math/Functions.h"
#include "Base/Util/Assert.h"
#include "Resample/Coherence/CoheringSubparticles.h"
#include "Resample/Element/DiffuseElement.h"
#include "Sample/Aggregate/InterferenceNone.h"
DecouplingApproximationStrategy::DecouplingApproximationStrategy(
const OwningVector<const CoheringSubparticles>& weighted_formfactors, const IInterference* iff,
SimulationOptions options, bool polarized)
: IInterparticleStrategy(weighted_formfactors, options, polarized)
, m_iff(iff ? iff->clone() : new InterferenceNone())
{
}
double DecouplingApproximationStrategy::scalarCalculation(const DiffuseElement& ele) const
{
double intensity = 0.0;
complex_t amplitude = 0;
for (const auto& ffw : m_weighted_formfactors) {
const complex_t ff = ffw->summedFF(ele);
ASSERT(std::isfinite(ff.real())); // numerical error in coherent sum?
const double fraction = ffw->relativeAbundance();
amplitude += fraction * ff;
intensity += fraction * std::norm(ff);
}
// Renaud, Lazzari, Leroy 2009, eq. (160).
// Incoherent cross section is intensity - norm(amplitude).
// Coherent cross section is S_q*norm(amplitude).
return intensity + (m_iff->structureFactor(ele.meanQ()) - 1) * std::norm(amplitude);
}
double DecouplingApproximationStrategy::polarizedCalculation(const DiffuseElement& ele) const
{
SpinMatrix mean_intensity;
SpinMatrix mean_amplitude;
const SpinMatrix& polarizer = ele.polarizer();
const SpinMatrix& analyzer = ele.analyzer();
for (const auto& ffw : m_weighted_formfactors) {
const SpinMatrix ff = ffw->summedPolFF(ele);
ASSERT(ff.allFinite()); // numerical error in coherent sum?
const double fraction = ffw->relativeAbundance();
mean_amplitude += fraction * ff;
mean_intensity += fraction * (ff * polarizer * ff.adjoint());
}
const SpinMatrix amplitude_matrix =
analyzer * mean_amplitude * polarizer * mean_amplitude.adjoint();
const SpinMatrix intensity_matrix = analyzer * mean_intensity;
const double amplitude_trace = std::abs(amplitude_matrix.trace());
const double intensity_trace = std::abs(intensity_matrix.trace());
return intensity_trace + (m_iff->structureFactor(ele.meanQ()) - 1) * amplitude_trace;
}
|
