diff --git a/Core/Algorithms/ScalarSpecularInfoMap.cpp b/Core/Algorithms/ScalarSpecularInfoMap.cpp
index 87a2cadcb0fad47a063746e82f45d2bca09edf1e..41dddfb610d443c85a76c338b3c257ada3b071e1 100644
--- a/Core/Algorithms/ScalarSpecularInfoMap.cpp
+++ b/Core/Algorithms/ScalarSpecularInfoMap.cpp
@@ -28,26 +28,26 @@ ScalarSpecularInfoMap *ScalarSpecularInfoMap::clone() const
return new ScalarSpecularInfoMap(mp_multilayer, m_layer);
}
-//! \todo Can we avoid the code duplication in the two following functions ?
-
const ScalarRTCoefficients *ScalarSpecularInfoMap::getOutCoefficients(
double alpha_f, double phi_f, double wavelength) const
{
- (void)phi_f;
- SpecularMatrix::MultiLayerCoeff_t coeffs;
- // phi has no effect on R,T, so just pass zero:
+ (void)phi_f; // phi has no effect on R,T, so just pass zero:
kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, -alpha_f, 0.0);
- SpecularMatrix::execute(*mp_multilayer, kvec, coeffs);
- return new ScalarRTCoefficients(coeffs[m_layer]);
+ return getCoefficients(kvec);
}
const ScalarRTCoefficients *ScalarSpecularInfoMap::getInCoefficients(
double alpha_i, double phi_i, double wavelength) const
{
- (void)phi_i;
- SpecularMatrix::MultiLayerCoeff_t coeffs;
- // phi has no effect on R,T, so just pass zero:
+ (void)phi_i; // phi has no effect on R,T, so just pass zero:
kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_i, 0.0);
+ return getCoefficients(kvec);
+}
+
+const ScalarRTCoefficients *ScalarSpecularInfoMap::getCoefficients(kvector_t kvec) const
+{
+ SpecularMatrix::MultiLayerCoeff_t coeffs;
SpecularMatrix::execute(*mp_multilayer, kvec, coeffs);
- return new ScalarRTCoefficients(coeffs[m_layer]);
+ auto layer_coeffs = coeffs[m_layer];
+ return new ScalarRTCoefficients(layer_coeffs);
}
diff --git a/Core/Algorithms/ScalarSpecularInfoMap.h b/Core/Algorithms/ScalarSpecularInfoMap.h
index deb364aafde344f512de63d67b54b147a6938dff..90cbaeb20220c60d18a1dc8194162a0aaed25dbb 100644
--- a/Core/Algorithms/ScalarSpecularInfoMap.h
+++ b/Core/Algorithms/ScalarSpecularInfoMap.h
@@ -45,6 +45,7 @@ public:
private:
const MultiLayer *mp_multilayer;
const int m_layer;
+ const ScalarRTCoefficients *getCoefficients(kvector_t kvec) const;
};
#endif /* SCALARSPECULARINFOMAP_H_ */
diff --git a/Core/Algorithms/SpecularMatrix.cpp b/Core/Algorithms/SpecularMatrix.cpp
index d158bb208a976f5ce560f0641485dcf0e91d3f4b..d6b33d70320ac857d71554e1a1d57e39ad484e01 100644
--- a/Core/Algorithms/SpecularMatrix.cpp
+++ b/Core/Algorithms/SpecularMatrix.cpp
@@ -24,6 +24,8 @@ namespace {
const complex_t imag_unit = complex_t(0.0, 1.0);
}
+void setZeroBelow(SpecularMatrix::MultiLayerCoeff_t& coeff, size_t current_layer);
+
// Computes refraction angles and transmission/reflection coefficients
// for given coherent wave propagation in a multilayer.
// k : length: wavenumber in vacuum, direction: defined in layer 0
@@ -59,9 +61,7 @@ void SpecularMatrix::execute(const MultiLayer& sample, const kvector_t k,
if (coeff[0].lambda == 0.0) {
coeff[0].t_r(0) = 1.0;
coeff[0].t_r(1) = -1.0;
- for (size_t i=1; i<N; ++i) {
- coeff[i].t_r.setZero();
- }
+ setZeroBelow(coeff, 0);
return;
}
// Calculate transmission/refraction coefficients t_r for each layer,
@@ -92,10 +92,32 @@ void SpecularMatrix::execute(const MultiLayer& sample, const kvector_t k,
(lambda_rough-lambda_below)*coeff[i+1].t_r(0) +
(lambda_rough+lambda_below)*coeff[i+1].t_r(1) )/2.0/lambda *
std::exp( ikd*lambda);
+ if (std::isinf(std::abs(coeff[i].getScalarT()))) {
+ coeff[i].t_r(0) = 1.0;
+ coeff[i].t_r(1) = 0.0;
+ setZeroBelow(coeff, i);
+ }
}
// Normalize to incoming downward travelling amplitude = 1.
complex_t T0 = coeff[0].getScalarT();
+ // This trick is used to avoid overflows in the intermediate steps of
+ // complex division:
+ double tmax = std::max(std::abs(T0.real()), std::abs(T0.imag()));
+ for (size_t i=0; i<N; ++i) {
+ coeff[i].t_r(0) /= tmax;
+ coeff[i].t_r(1) /= tmax;
+ }
+ // Now the real normalization to T_0 proceeds
+ T0 = coeff[0].getScalarT();
for (size_t i=0; i<N; ++i) {
coeff[i].t_r = coeff[i].t_r/T0;
}
}
+
+void setZeroBelow(SpecularMatrix::MultiLayerCoeff_t& coeff, size_t current_layer)
+{
+ size_t N = coeff.size();
+ for (size_t i=current_layer+1; i<N; ++i) {
+ coeff[i].t_r.setZero();
+ }
+}
diff --git a/Core/FormFactors/FormFactorDWBA.cpp b/Core/FormFactors/FormFactorDWBA.cpp
index c0b389bf64bb5cde187e9d37070285a7d2df3ed7..7a829b95a87025b798e967dce61361e2fe883c43 100644
--- a/Core/FormFactors/FormFactorDWBA.cpp
+++ b/Core/FormFactors/FormFactorDWBA.cpp
@@ -47,7 +47,8 @@ void FormFactorDWBA::setSpecularInfo(const ILayerRTCoefficients *p_in_coeffs,
complex_t FormFactorDWBA::evaluate(const WavevectorInfo& wavevectors) const
{
calculateTerms(wavevectors);
- return m_term_S + m_term_RS + m_term_SR + m_term_RSR;
+ complex_t sum_terms = m_term_S + m_term_RS + m_term_SR + m_term_RSR;
+ return sum_terms;
}
void FormFactorDWBA::calculateTerms(const WavevectorInfo& wavevectors) const
@@ -71,14 +72,16 @@ void FormFactorDWBA::calculateTerms(const WavevectorInfo& wavevectors) const
WavevectorInfo k_TR(k_i_T, k_f_R, wavelength);
WavevectorInfo k_RR(k_i_R, k_f_R, wavelength);
+ // Get the four R,T coefficients
+ complex_t T_in = mp_in_coeffs->getScalarT();
+ complex_t R_in = mp_in_coeffs->getScalarR();
+ complex_t T_out = mp_out_coeffs->getScalarT();
+ complex_t R_out = mp_out_coeffs->getScalarR();
+
// The four different scattering contributions; S stands for scattering
// off the particle, R for reflection off the layer interface
- m_term_S = mp_in_coeffs->getScalarT()*mp_form_factor->evaluate(k_TT)
- * mp_out_coeffs->getScalarT();
- m_term_RS = mp_in_coeffs->getScalarR()*mp_form_factor->evaluate(k_RT)
- * mp_out_coeffs->getScalarT();
- m_term_SR = mp_in_coeffs->getScalarT()*mp_form_factor->evaluate(k_TR)
- * mp_out_coeffs->getScalarR();
- m_term_RSR = mp_in_coeffs->getScalarR()*mp_form_factor->evaluate(k_RR)
- * mp_out_coeffs->getScalarR();
+ m_term_S = T_in * mp_form_factor->evaluate(k_TT) * T_out;
+ m_term_RS = R_in * mp_form_factor->evaluate(k_RT) * T_out;
+ m_term_SR = T_in * mp_form_factor->evaluate(k_TR) * R_out;
+ m_term_RSR = R_in * mp_form_factor->evaluate(k_RR) * R_out;
}