From 0246b57c2cf76bdd5fdf683b3e830a200f6e02c1 Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Mon, 6 May 2013 14:00:30 +0200
Subject: [PATCH] Refactored SpecularMatrix to use Eigen library for matrix
calculations
---
Core/Algorithms/inc/SpecularMatrix.h | 23 ++++++--------
Core/Algorithms/src/SpecularMatrix.cpp | 43 ++++++++------------------
shared.pri | 7 +++++
3 files changed, 29 insertions(+), 44 deletions(-)
diff --git a/Core/Algorithms/inc/SpecularMatrix.h b/Core/Algorithms/inc/SpecularMatrix.h
index 7913bf67867..70b8cec3580 100644
--- a/Core/Algorithms/inc/SpecularMatrix.h
+++ b/Core/Algorithms/inc/SpecularMatrix.h
@@ -20,29 +20,27 @@
#include "ISimulation.h"
#include "MultiLayer.h"
+#include "Eigen/Core"
+
//! Implements the matrix formalism for the calculation of wave amplitudes of
//! the coherent wave solution in a multilayer
class SpecularMatrix : public ISimulation
{
public:
- SpecularMatrix() : m_use_roughness(false) {}
+ SpecularMatrix() : m_use_roughness(false) { (void)m_use_roughness; }
//! layer coefficients for matrix formalism
class LayerMatrixCoeff {
public:
- LayerMatrixCoeff() : lambda(0), phi(0), psi(0), l11(0), l12(0), l21(0), l22(0) {}
+ LayerMatrixCoeff() : lambda(0) {}
~LayerMatrixCoeff() {}
- complex_t R() const { return (psi+phi/lambda)/2.0; }
- complex_t T() const { return (psi-phi/lambda)/2.0; }
+ complex_t R() const { return (phi_psi(1)+phi_psi(0)/lambda)/2.0; }
+ complex_t T() const { return (phi_psi(1)-phi_psi(0)/lambda)/2.0; }
// A - amplitude of initial wave, R, T - amplitudes of reflected and transmitted waves
complex_t lambda; // positive eigenvalue of transfer matrix
- complex_t phi; // amplitude of the normalized derivative of wavefunction at top
- complex_t psi; // amplitude of the wavefunction at top
- complex_t l11; // matrix element of transfer matrix
- complex_t l12; // matrix element of transfer matrix
- complex_t l21; // matrix element of transfer matrix
- complex_t l22; // matrix element of transfer matrix
+ Eigen::Vector2cd phi_psi;
+ Eigen::Matrix2cd l;
};
//! multi layer coefficients for matrix formalism
@@ -55,10 +53,7 @@ public:
inline void clear() { m_data.clear(); }
inline void resize(size_t size) { m_data.resize(size); }
complex_t R; // total reflection coefficient
- complex_t L11;
- complex_t L12;
- complex_t L21;
- complex_t L22;
+ Eigen::Matrix2cd L;
private:
std::vector<LayerMatrixCoeff > m_data;
};
diff --git a/Core/Algorithms/src/SpecularMatrix.cpp b/Core/Algorithms/src/SpecularMatrix.cpp
index 28ca5ce0aff..f078ad3790c 100644
--- a/Core/Algorithms/src/SpecularMatrix.cpp
+++ b/Core/Algorithms/src/SpecularMatrix.cpp
@@ -47,37 +47,21 @@ void SpecularMatrix::calculateTransferMatrices(const MultiLayer& sample,
{
double ksinalpha = k.mag()*std::abs( k.cosTheta() );
// Layer 0 gets identity matrix:
- coeff[0].l11 = complex_t(1.0, 0.0);
- coeff[0].l12 = complex_t(0.0, 0.0);
- coeff[0].l21 = complex_t(0.0, 0.0);
- coeff[0].l22 = complex_t(1.0, 0.0);
- coeff.L11 = complex_t(1.0, 0.0);
- coeff.L12 = complex_t(0.0, 0.0);
- coeff.L21 = complex_t(0.0, 0.0);
- coeff.L22 = complex_t(1.0, 0.0);
+ coeff[0].l.setIdentity();
+ coeff.L.setIdentity();
for(size_t i=1; i<coeff.size()-1; ++i) {
complex_t lambda = coeff[i].lambda;
complex_t coskdsinlambda = std::cos(ksinalpha*sample.getLayer(i)->getThickness()*lambda);
complex_t sinkdsinlambda = std::sin(ksinalpha*sample.getLayer(i)->getThickness()*lambda);
- coeff[i].l11 = coskdsinlambda;
- coeff[i].l12 = -complex_t(0.0, 1.0)*lambda*sinkdsinlambda;
- coeff[i].l21 = -complex_t(0.0, 1.0)*sinkdsinlambda/lambda;
- coeff[i].l22 = coskdsinlambda;
- complex_t ltemp11 = coeff[i].l11*coeff.L11 + coeff[i].l12*coeff.L21;
- complex_t ltemp12 = coeff[i].l11*coeff.L12 + coeff[i].l12*coeff.L22;
- complex_t ltemp21 = coeff[i].l21*coeff.L11 + coeff[i].l22*coeff.L21;
- complex_t ltemp22 = coeff[i].l21*coeff.L12 + coeff[i].l22*coeff.L22;
- coeff.L11 = ltemp11;
- coeff.L12 = ltemp12;
- coeff.L21 = ltemp21;
- coeff.L22 = ltemp22;
+ coeff[i].l(0,0) = coskdsinlambda;
+ coeff[i].l(0,1) = -complex_t(0.0, 1.0)*lambda*sinkdsinlambda;
+ coeff[i].l(1,0) = -complex_t(0.0, 1.0)*sinkdsinlambda/lambda;
+ coeff[i].l(1,1) = coskdsinlambda;
+ coeff.L = coeff[i].l * coeff.L;
}
// Last layer also gets identity matrix:
size_t N = coeff.size();
- coeff[N-1].l11 = complex_t(1.0, 0.0);
- coeff[N-1].l12 = complex_t(0.0, 0.0);
- coeff[N-1].l21 = complex_t(0.0, 0.0);
- coeff[N-1].l22 = complex_t(1.0, 0.0);
+ coeff[N-1].l.setIdentity();
}
void SpecularMatrix::calculateBoundaryValues(MultiLayerCoeff_t& coeff) const
@@ -85,14 +69,13 @@ void SpecularMatrix::calculateBoundaryValues(MultiLayerCoeff_t& coeff) const
complex_t lambda0 = coeff[0].lambda;
size_t N = coeff.size();
complex_t lambdaN = coeff[N-1].lambda;
- complex_t denominator = (lambda0*coeff.L11 + coeff.L12 + lambdaN*(lambda0*coeff.L21 + coeff.L22));
- coeff.R = (lambda0*coeff.L11 - coeff.L12 + lambdaN*(lambda0*coeff.L21 - coeff.L22))/denominator;
+ complex_t denominator = (lambda0*coeff.L(0,0) + coeff.L(0,1) + lambdaN*(lambda0*coeff.L(1,0) + coeff.L(1,1)));
+ coeff.R = (lambda0*coeff.L(0,0) - coeff.L(0,1) + lambdaN*(lambda0*coeff.L(1,0) - coeff.L(1,1)))/denominator;
// Boundary values at bottom of top layer
// and top of first layer:
- coeff[1].phi = coeff[0].phi = lambda0*(coeff.R - 1.0);
- coeff[1].psi = coeff[0].psi = coeff.R + 1.0;
+ coeff[1].phi_psi(0) = coeff[0].phi_psi(0) = lambda0*(coeff.R - 1.0);
+ coeff[1].phi_psi(1) = coeff[0].phi_psi(1) = coeff.R + 1.0;
for(size_t i=2; i<coeff.size(); ++i) {
- coeff[i].phi = coeff[i-1].l11*coeff[i-1].phi + coeff[i-1].l12*coeff[i-1].psi;
- coeff[i].psi = coeff[i-1].l21*coeff[i-1].phi + coeff[i-1].l22*coeff[i-1].psi;
+ coeff[i].phi_psi = coeff[i-1].l*coeff[i-1].phi_psi;
}
}
diff --git a/shared.pri b/shared.pri
index bfa23a29031..cfa5c631b8e 100644
--- a/shared.pri
+++ b/shared.pri
@@ -33,6 +33,13 @@ INCLUDEPATH *= $${GSL_INCLUDE}
LIBS *= -L$${GSL_LIB}
LIBS += -lgsl -lgslcblas
+# --- checking eigen headers ---
+EIGEN_HEADERFILE = Eigen/Core
+EIGEN_HEADER_LOCATIONS = /opt/local/include /usr/local/include /usr/include
+for(dir, EIGEN_HEADER_LOCATIONS): isEmpty(EIGEN_INCLUDE): exists($${dir}/$${EIGEN_HEADERFILE}): EIGEN_INCLUDE = $${dir}
+isEmpty(EIGEN_INCLUDE): message("Can't find" $${EIGEN_HEADERFILE} "in" $${EIGEN_HEADER_LOCATIONS})
+INCLUDEPATH *= $${EIGEN_INCLUDE}
+
# --- checking fftw3 ---
FFTW3_HEADERFILE = fftw3.h
FFTW3_HEADER_LOCATIONS = /opt/local/include /usr/local/include /usr/include
--
GitLab