From b889d91c4dcdfb9370c72a7aa28dd010daa9f42b Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (h)" <j.wuttke@fz-juelich.de>
Date: Sun, 21 Feb 2016 14:56:29 +0100
Subject: [PATCH] Prompted by needs in polyhedral FF computation, rework
 BasisVector3D. In particular, member function setLambdaAlphaPhi -> static
 vecOfLambdaAlphaPhi.

---
 Core/Algorithms/inc/LayerDWBASimulation.h     |   3 +-
 Core/Algorithms/src/Beam.cpp                  |   4 +-
 .../src/IInterferenceFunctionStrategy.cpp     |  38 ++--
 Core/Algorithms/src/LayerDWBASimulation.cpp   |   4 +-
 Core/Algorithms/src/MatrixSpecularInfoMap.cpp |   6 +-
 Core/Algorithms/src/RectangularDetector.cpp   |  10 +-
 Core/Algorithms/src/ScalarSpecularInfoMap.cpp |   8 +-
 Core/Algorithms/src/SimulationElement.cpp     |   4 +-
 Core/Algorithms/src/SpecularSimulation.cpp    |   3 +-
 Core/Algorithms/src/SphericalDetector.cpp     |   4 +-
 Core/Geometry/inc/BasicVector3D.h             | 173 ++++++++----------
 Core/Geometry/src/BasicVector3D.cpp           |  21 +++
 Core/Samples/src/Lattice.cpp                  |  14 +-
 Core/Tools/inc/Bin.h                          |  12 +-
 14 files changed, 139 insertions(+), 165 deletions(-)

diff --git a/Core/Algorithms/inc/LayerDWBASimulation.h b/Core/Algorithms/inc/LayerDWBASimulation.h
index f1ea2ab043f..8db18acee60 100644
--- a/Core/Algorithms/inc/LayerDWBASimulation.h
+++ b/Core/Algorithms/inc/LayerDWBASimulation.h
@@ -38,8 +38,7 @@ public:
     void setSpecularInfo(const LayerSpecularInfo &specular_info);
 
 protected:
-    Bin1DCVector getKfBin(double wavelength, const Bin1D& alpha_bin,
-                          const Bin1D& phi_bin) const;
+    Bin1DCVector getKfBin(double wavelength, const Bin1D& alpha_bin, const Bin1D& phi_bin) const;
     Layer *mp_layer;
     LayerSpecularInfo *mp_specular_info;
 };
diff --git a/Core/Algorithms/src/Beam.cpp b/Core/Algorithms/src/Beam.cpp
index 9cb792b63bb..59fdfa24ada 100644
--- a/Core/Algorithms/src/Beam.cpp
+++ b/Core/Algorithms/src/Beam.cpp
@@ -47,9 +47,7 @@ Beam &Beam::operator=(const Beam &other)
 
 kvector_t Beam::getCentralK() const
 {
-    kvector_t k;
-    k.setLambdaAlphaPhi(m_wavelength, -1.0 * m_alpha, m_phi);
-    return k;
+    return Geometry::vecOfLambdaAlphaPhi(m_wavelength, -m_alpha, m_phi);
 }
 
 void Beam::setCentralK(double wavelength, double alpha_i, double phi_i)
diff --git a/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp b/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp
index ab28d28dd47..cd284503d9c 100644
--- a/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp
+++ b/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp
@@ -27,10 +27,6 @@ IInterferenceFunctionStrategy::IInterferenceFunctionStrategy(SimulationParameter
                 this, &IInterferenceFunctionStrategy::evaluate_for_fixed_angles_pol, 2);
 }
 
-IInterferenceFunctionStrategy::~IInterferenceFunctionStrategy()
-{
-}
-
 void IInterferenceFunctionStrategy::init(const SafePointerVector<FormFactorInfo> &form_factor_infos,
                                          const IInterferenceFunction& iff)
 {
@@ -70,23 +66,20 @@ void IInterferenceFunctionStrategy::calculateFormFactorList(
     clearFormFactorLists();
 
     double wavelength = sim_element.getWavelength();
-    complex_t wavevector_scattering_factor = Units::PI/wavelength/wavelength;
+    double wavevector_scattering_factor = Units::PI/wavelength/wavelength;
     double alpha_i = sim_element.getAlphaI();
     double phi_i = sim_element.getPhiI();
-    cvector_t k_i;
-    k_i.setLambdaAlphaPhi(wavelength, alpha_i, phi_i);
+    cvector_t k_i = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_i, phi_i).complex();
     WavevectorInfo wavevectors(k_i, Geometry::toComplexVector(sim_element.getMeanKF()), wavelength);
 
     boost::scoped_ptr<const ILayerRTCoefficients> P_in_coeffs(
         mP_specular_info->getInCoefficients(alpha_i, 0.0, wavelength));
     boost::scoped_ptr<const ILayerRTCoefficients> P_out_coeffs(
         mP_specular_info->getOutCoefficients(sim_element.getAlphaMean(), 0.0, wavelength));
-    SafePointerVector<FormFactorInfo>::const_iterator it = m_ff_infos.begin();
-    while (it != m_ff_infos.end()) {
-        (*it)->mp_ff->setSpecularInfo(P_in_coeffs.get(), P_out_coeffs.get());
-        complex_t ff_mat = (*it)->mp_ff->evaluate(wavevectors);
+    for( auto it: m_ff_infos ) {
+        it->mp_ff->setSpecularInfo(P_in_coeffs.get(), P_out_coeffs.get());
+        complex_t ff_mat = it->mp_ff->evaluate(wavevectors);
         m_ff.push_back(wavevector_scattering_factor*ff_mat);
-        ++it;
     }
 }
 
@@ -96,11 +89,10 @@ void IInterferenceFunctionStrategy::calculateFormFactorLists(
     clearFormFactorLists();
 
     double wavelength = sim_element.getWavelength();
-    complex_t wavevector_scattering_factor = Units::PI/wavelength/wavelength;
+    double wavevector_scattering_factor = Units::PI/wavelength/wavelength;
     double alpha_i = sim_element.getAlphaI();
     double phi_i = sim_element.getPhiI();
-    cvector_t k_i;
-    k_i.setLambdaAlphaPhi(wavelength, alpha_i, phi_i);
+    cvector_t k_i = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_i, phi_i).complex();
     WavevectorInfo wavevectors(k_i, Geometry::toComplexVector(sim_element.getMeanKF()), wavelength);
 
     boost::scoped_ptr<const ILayerRTCoefficients> P_in_coeffs(
@@ -108,12 +100,10 @@ void IInterferenceFunctionStrategy::calculateFormFactorLists(
     boost::scoped_ptr<const ILayerRTCoefficients> P_out_coeffs(
         mP_specular_info->getOutCoefficients(sim_element.getAlphaMean(), sim_element.getPhiMean(),
                                              wavelength));
-    SafePointerVector<FormFactorInfo>::const_iterator it = m_ff_infos.begin();
-    while (it != m_ff_infos.end()) {
-        (*it)->mp_ff->setSpecularInfo(P_in_coeffs.get(), P_out_coeffs.get());
-        Eigen::Matrix2cd ff_mat = (*it)->mp_ff->evaluatePol(wavevectors);
+    for ( auto it: m_ff_infos ) {
+        it->mp_ff->setSpecularInfo(P_in_coeffs.get(), P_out_coeffs.get());
+        Eigen::Matrix2cd ff_mat = it->mp_ff->evaluatePol(wavevectors);
         m_ff_pol.push_back(wavevector_scattering_factor*ff_mat);
-        ++it;
     }
 }
 
@@ -127,8 +117,8 @@ double IInterferenceFunctionStrategy::MCIntegratedEvaluate(const SimulationEleme
 {
     double min_array[] = {0.0, 0.0};
     double max_array[] = {1.0, 1.0};
-    return mP_integrator->integrate(min_array, max_array, (void *)&sim_element,
-                                            m_sim_params.m_mc_points);
+    return mP_integrator->integrate(
+        min_array, max_array, (void *)&sim_element, m_sim_params.m_mc_points);
 }
 
 double IInterferenceFunctionStrategy::MCIntegratedEvaluatePol(
@@ -136,8 +126,8 @@ double IInterferenceFunctionStrategy::MCIntegratedEvaluatePol(
 {
     double min_array[] = {0.0, 0.0};
     double max_array[] = {1.0, 1.0};
-    return mP_integrator_pol->integrate(min_array, max_array, (void *)&sim_element,
-                                            m_sim_params.m_mc_points);
+    return mP_integrator_pol->integrate(
+        min_array, max_array, (void *)&sim_element, m_sim_params.m_mc_points);
 }
 
 double IInterferenceFunctionStrategy::evaluate_for_fixed_angles(double *fractions, size_t dim,
diff --git a/Core/Algorithms/src/LayerDWBASimulation.cpp b/Core/Algorithms/src/LayerDWBASimulation.cpp
index d78f4eaf692..8725a83a9d5 100644
--- a/Core/Algorithms/src/LayerDWBASimulation.cpp
+++ b/Core/Algorithms/src/LayerDWBASimulation.cpp
@@ -30,8 +30,8 @@ LayerDWBASimulation::~LayerDWBASimulation()
     delete mp_specular_info;
 }
 
-Bin1DCVector LayerDWBASimulation::getKfBin(double wavelength,
-        const Bin1D& alpha_bin, const Bin1D& phi_bin) const
+Bin1DCVector LayerDWBASimulation::getKfBin(
+    double wavelength, const Bin1D& alpha_bin, const Bin1D& phi_bin) const
 {
     Bin1DCVector k_f_bin(wavelength, alpha_bin, phi_bin);
     return k_f_bin;
diff --git a/Core/Algorithms/src/MatrixSpecularInfoMap.cpp b/Core/Algorithms/src/MatrixSpecularInfoMap.cpp
index 7b81a8e1215..feea5b6cb29 100644
--- a/Core/Algorithms/src/MatrixSpecularInfoMap.cpp
+++ b/Core/Algorithms/src/MatrixSpecularInfoMap.cpp
@@ -39,8 +39,7 @@ const MatrixRTCoefficients *MatrixSpecularInfoMap::getOutCoefficients(
         double alpha_f, double phi_f, double wavelength) const
 {
     SpecularMagnetic::MultiLayerCoeff_t coeffs;
-    kvector_t kvec;
-    kvec.setLambdaAlphaPhi(wavelength, alpha_f, phi_f);
+    kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_f, phi_f);
     SpecularMagnetic::execute(*mP_inverted_multilayer, -kvec, coeffs);
     return new MatrixRTCoefficients(coeffs[m_layer]);
 }
@@ -49,8 +48,7 @@ const MatrixRTCoefficients *MatrixSpecularInfoMap::getInCoefficients(
         double alpha_i, double phi_i, double wavelength) const
 {
     SpecularMagnetic::MultiLayerCoeff_t coeffs;
-    kvector_t kvec;
-    kvec.setLambdaAlphaPhi(wavelength, alpha_i, phi_i);
+    kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_i, phi_i);
     SpecularMagnetic::execute(*mP_multilayer, kvec, coeffs);
     return new MatrixRTCoefficients(coeffs[m_layer]);
 }
diff --git a/Core/Algorithms/src/RectangularDetector.cpp b/Core/Algorithms/src/RectangularDetector.cpp
index 6a98fe78463..5b4b49a2765 100644
--- a/Core/Algorithms/src/RectangularDetector.cpp
+++ b/Core/Algorithms/src/RectangularDetector.cpp
@@ -132,9 +132,8 @@ IPixelMap *RectangularDetector::createPixelMap(size_t index) const
     return new RectPixelMap(corner_position, width, height);
 }
 
-std::string RectangularDetector::addParametersToExternalPool(std::string path,
-                                                             ParameterPool *external_pool,
-                                                             int copy_number) const
+std::string RectangularDetector::addParametersToExternalPool(
+    std::string path, ParameterPool *external_pool, int copy_number) const
 {
     // add own parameters
     std::string new_path
@@ -275,7 +274,8 @@ OutputData<double> *RectangularDetector::createDetectorMap(const Beam &beam,
 std::vector<IDetector2D::EAxesUnits> RectangularDetector::getValidAxesUnits() const
 {
     std::vector<IDetector2D::EAxesUnits> result = IDetector2D::getValidAxesUnits();
-    std::vector<IDetector2D::EAxesUnits> addon = {IDetector2D::RADIANS, IDetector2D::DEGREES, IDetector2D::MM, IDetector2D::QYQZ};
+    std::vector<IDetector2D::EAxesUnits> addon =
+        { IDetector2D::RADIANS, IDetector2D::DEGREES, IDetector2D::MM, IDetector2D::QYQZ };
     result.insert(result.end(), addon.begin(), addon.end());
     return result;
 }
@@ -365,13 +365,11 @@ void RectangularDetector::initNormalVector(const kvector_t &central_k)
     else if(m_detector_arrangement == PERPENDICULAR_TO_REFLECTED_BEAM_DPOS) {
         m_normal_to_detector = m_distance*central_k_unit;
         m_normal_to_detector.setZ(-m_normal_to_detector.z());
-
     }
 
     else {
         throw LogicErrorException("RectangularDetector::init() -> Unknown detector arrangement");
     }
-
 }
 
 void RectangularDetector::initUandV(double alpha_i)
diff --git a/Core/Algorithms/src/ScalarSpecularInfoMap.cpp b/Core/Algorithms/src/ScalarSpecularInfoMap.cpp
index a7b0aa178d1..5e4f296686f 100644
--- a/Core/Algorithms/src/ScalarSpecularInfoMap.cpp
+++ b/Core/Algorithms/src/ScalarSpecularInfoMap.cpp
@@ -28,14 +28,15 @@ 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;
-    kvector_t kvec;
     // phi has no effect on R,T, so just pass zero:
-    kvec.setLambdaAlphaPhi(wavelength, -alpha_f, 0.0);
+    kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, -alpha_f, 0.0);
     SpecularMatrix::execute(*mp_multilayer, kvec, coeffs);
     return new ScalarRTCoefficients(coeffs[m_layer]);
 }
@@ -45,9 +46,8 @@ const ScalarRTCoefficients *ScalarSpecularInfoMap::getInCoefficients(
 {
     (void)phi_i;
     SpecularMatrix::MultiLayerCoeff_t coeffs;
-    kvector_t kvec;
     // phi has no effect on R,T, so just pass zero:
-    kvec.setLambdaAlphaPhi(wavelength, alpha_i, 0.0);
+    kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_i, 0.0);
     SpecularMatrix::execute(*mp_multilayer, kvec, coeffs);
     return new ScalarRTCoefficients(coeffs[m_layer]);
 }
diff --git a/Core/Algorithms/src/SimulationElement.cpp b/Core/Algorithms/src/SimulationElement.cpp
index 79cbe096ebd..50c564b873d 100644
--- a/Core/Algorithms/src/SimulationElement.cpp
+++ b/Core/Algorithms/src/SimulationElement.cpp
@@ -55,9 +55,7 @@ SimulationElement::SimulationElement(const SimulationElement &other, double x, d
 
 kvector_t SimulationElement::getKI() const
 {
-    kvector_t k_i;
-    k_i.setLambdaAlphaPhi(m_wavelength, m_alpha_i, m_phi_i);
-    return k_i;
+    return Geometry::vecOfLambdaAlphaPhi(m_wavelength, m_alpha_i, m_phi_i);
 }
 
 kvector_t SimulationElement::getMeanKF() const
diff --git a/Core/Algorithms/src/SpecularSimulation.cpp b/Core/Algorithms/src/SpecularSimulation.cpp
index ba82bf64086..c9188d3b8c9 100644
--- a/Core/Algorithms/src/SpecularSimulation.cpp
+++ b/Core/Algorithms/src/SpecularSimulation.cpp
@@ -242,8 +242,7 @@ void SpecularSimulation::collectRTCoefficientsScalar(const MultiLayer *multilaye
     OutputData<MultiLayerRTCoefficients_t>::iterator it = m_data.begin();
     while (it != m_data.end()) {
         double alpha_i = m_data.getAxisValue(it.getIndex(), 0);
-        kvector_t kvec;
-        kvec.setLambdaAlphaPhi(m_lambda, -alpha_i, 0.0);
+        kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(m_lambda, -alpha_i, 0.0);
 
         SpecularMatrix::MultiLayerCoeff_t coeffs;
         SpecularMatrix::execute(*multilayer, kvec, coeffs);
diff --git a/Core/Algorithms/src/SphericalDetector.cpp b/Core/Algorithms/src/SphericalDetector.cpp
index 6e77f8f9cdd..0e50dac486f 100644
--- a/Core/Algorithms/src/SphericalDetector.cpp
+++ b/Core/Algorithms/src/SphericalDetector.cpp
@@ -224,11 +224,9 @@ AngularPixelMap *AngularPixelMap::createZeroSizeMap(double x, double y) const
 
 kvector_t AngularPixelMap::getK(double x, double y, double wavelength) const
 {
-    kvector_t result;
     double phi = m_phi + x*m_dphi;
     double alpha = m_alpha + y*m_dalpha;
-    result.setLambdaAlphaPhi(wavelength, alpha, phi);
-    return result;
+    return Geometry::vecOfLambdaAlphaPhi(wavelength, alpha, phi);
 }
 
 double AngularPixelMap::getIntegrationFactor(double x, double y) const
diff --git a/Core/Geometry/inc/BasicVector3D.h b/Core/Geometry/inc/BasicVector3D.h
index 4ddb49ff3ed..e64735cfb6d 100644
--- a/Core/Geometry/inc/BasicVector3D.h
+++ b/Core/Geometry/inc/BasicVector3D.h
@@ -6,9 +6,8 @@
 //! @brief      Defines template class BasicVector3D.
 //!
 //! Forked from CLHEP/Geometry by E. Chernyaev <Evgueni.Tcherniaev@cern.ch>,
-//! then reduced to rotations and mostly rewritten; point and vector semantics
-//! is no longer represented by class type; transforms are no longer methods of
-//! BasicVector3D; there is a new class Transform3D.
+//! then reworked beyond recongnition. Removed split of point and vector semantics.
+//! Transforms are relegated to a separate class Transform3D.
 //!
 //! @homepage  http://bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
@@ -29,12 +28,7 @@ namespace Geometry {
 
 //! @class BasicVector3D
 //! @ingroup tools_internal
-//! @brief Base class for Point3D<T>, Vector3D<T> and Normal3D<T>.
-//!
-//! It defines only common functionality for those classes and
-//! should not be used as separate class.
-//!
-//! @author Evgeni Chernyaev 1996-2003
+//! @brief Three-dimensional vector template, for use with integer, double, or complex components.
 
 template<class T>
 class BasicVector3D {
@@ -42,16 +36,20 @@ protected:
     T v_[3];
 
 public:
+    // -------------------------------------------------------------------------
+    // Constructors and other set functions
+    // -------------------------------------------------------------------------
+
     //! Default constructor.
     BasicVector3D()
     { v_[0] = 0.0; v_[1] = 0.0; v_[2] = 0.0; }
 
-    //! Constructor from three numbers.
+    //! Constructor from cartesian components.
     BasicVector3D(const T x1, const T y1, const T z1)
     { v_[0] = x1; v_[1] = y1; v_[2] = z1; }
 
     // -------------------------------------------------------------------------
-    // Read and set components
+    // Component access
     // -------------------------------------------------------------------------
 
     //! Returns components by index.
@@ -99,7 +97,7 @@ public:
     { v_[0] /= a; v_[1] /= a; v_[2] /= a; return *this; }
 
     // -------------------------------------------------------------------------
-    // Norm
+    // Functions of this (with no further argument)
     // -------------------------------------------------------------------------
 
     //! Returns squared magnitude squared of the vector.
@@ -108,10 +106,6 @@ public:
     //! Returns magnitude of the vector.
     T mag() const;
 
-    // -------------------------------------------------------------------------
-    // Cylindrical and spherical coordinate systems
-    // -------------------------------------------------------------------------
-
     //! Returns squared distance from z axis.
     T magxy2() const;
 
@@ -130,8 +124,19 @@ public:
     //! Returns squared sine of polar angle.
     double sin2Theta() const;
 
+    //! Returns unit vector in direction of this (or null vector).
+    inline BasicVector3D<T> unit() const {
+        double len = mag();
+        return (len > 0.0) ?
+            BasicVector3D<T>(x()/len, y()/len, z()/len) :
+            BasicVector3D<T>();
+    }
+
+    //! Returns this, trivially converted to complex type.
+    BasicVector3D<std::complex<double>> complex() const;
+
     // -------------------------------------------------------------------------
-    // Combine two vectors
+    // Functions of this and another vector
     // -------------------------------------------------------------------------
 
     //! Returns dot product of vectors (antilinear in the first [=self] argument).
@@ -140,19 +145,22 @@ public:
     //! Returns cross product of vectors.
     BasicVector3D<T> cross(const BasicVector3D<T>& v ) const;
 
-    //! Returns normalized vector
-    BasicVector3D<T> normalize() const;
-
     //! Returns square of transverse component with respect to given axis.
-    double perp2(const BasicVector3D<T>& v) const;
+    //! Only for T=double.
+    double perp2(const BasicVector3D<double>& v) const;
 
     //! Returns transverse component with respect to given axis.
-    inline T perp(const BasicVector3D<T>& v) const
+    //! Only for T=double.
+    inline double perp(const BasicVector3D<double>& v) const
     { return std::sqrt(perp2(v)); }
 
     //! Returns angle with respect to another vector.
     double angle(const BasicVector3D<T>& v) const;
 
+    //! Returns projection of this onto other vector: (this*v)*v/|v|^2.
+    inline BasicVector3D<T> project(const BasicVector3D<T>& v) const
+    { return dot(v)*v/v.mag2(); }
+
     // -------------------------------------------------------------------------
     // Rotations
     // -------------------------------------------------------------------------
@@ -165,36 +173,11 @@ public:
     BasicVector3D<T> rotatedZ(double a) const;
     //! Returns result of rotation around the axis specified by another vector.
     BasicVector3D<T> rotated(double a, const BasicVector3D<T>& v) const;
-
-    // -------------------------------------------------------------------------
-    // Related vectors
-    // -------------------------------------------------------------------------
-
-    //! Returns unit vector in direction of this (or null vector).
-    inline BasicVector3D<T> unit() const {
-        double len = mag();
-        return (len > 0.0) ?
-            BasicVector3D<T>(x()/len, y()/len, z()/len) :
-            BasicVector3D<T>();
-    }
-
-    // -------------------------------------------------------------------------
-    // Specifically for grazing-incidence scattering
-    // -------------------------------------------------------------------------
-
-    //! Sets wave vector for given wavelength and angles/
-    inline void setLambdaAlphaPhi(
-        const T& _lambda, const T& _alpha, const T& _phi)
-        {
-            T k = PI2/_lambda;
-            v_[0] = k*std::cos(_alpha) * std::cos(_phi);
-            v_[1] = -k*std::cos(_alpha) * std::sin(_phi);
-            v_[2] = k*std::sin(_alpha);
-        }
 };
 
+
 // =============================================================================
-// Non-member functions for BasicVector3D<T>
+// Non-member functions
 // =============================================================================
 
 //! Output to stream.
@@ -205,94 +188,85 @@ operator<<(std::ostream& os, const BasicVector3D<T>& a)
 { return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")"; }
 
 // -----------------------------------------------------------------------------
-// Scalar operations
+// Unary operators
 // -----------------------------------------------------------------------------
 
 //! Unary plus.
 //! @relates BasicVector3D
 template <class T>
-inline BasicVector3D<T>
-operator+(const BasicVector3D<T>& v)
+inline BasicVector3D<T> operator+ (const BasicVector3D<T>& v)
 { return v; }
 
 //! Unary minus.
 //! @relates BasicVector3D
 template <class T>
-inline BasicVector3D<T>
-operator-(const BasicVector3D<T>& v)
+inline BasicVector3D<T> operator- (const BasicVector3D<T>& v)
 { return BasicVector3D<T>(-v.x(), -v.y(), -v.z()); }
 
-//! Multiplication vector by scalar.
-//! @relates BasicVector3D
-template <class T, class U>
-inline BasicVector3D<T>
-operator*(const BasicVector3D<T>& v, U a)
-{ return BasicVector3D<T>(v.x()*a, v.y()*a, v.z()*a); }
-
-//! Multiplication scalar by vector.
-//! @relates BasicVector3D
-template <class T, class U>
-inline BasicVector3D<T>
-operator*(U a, const BasicVector3D<T>& v)
-{ return BasicVector3D<T>(a*v.x(), a*v.y(), a*v.z()); }
-
-//! Division vector by scalar.
-//! @relates BasicVector3D
-template <class T, class U>
-inline BasicVector3D<T>
-operator/(const BasicVector3D<T>& v, U a)
-{ return BasicVector3D<T>(v.x()/a, v.y()/a, v.z()/a); }
-
-// -----------------------------------
-// Operations involving another vector
-// -----------------------------------
+// -----------------------------------------------------------------------------
+// Binary operators
+// -----------------------------------------------------------------------------
 
 //! Addition of two vectors.
 //! @relates BasicVector3D
 template <class T>
-inline BasicVector3D<T>
-operator+(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
+inline BasicVector3D<T> operator+(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
 { return BasicVector3D<T>(a.x()+b.x(), a.y()+b.y(), a.z()+b.z()); }
 
 //! Subtraction of two vectors.
 //! @relates BasicVector3D
 template <class T>
-inline BasicVector3D<T>
-operator-(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
+inline BasicVector3D<T> operator-(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
 { return BasicVector3D<T>(a.x()-b.x(), a.y()-b.y(), a.z()-b.z()); }
 
-//! Scalar product of two vectors.
+//! Multiplication vector by scalar.
+//! @relates BasicVector3D
+template <class T, class U>
+inline BasicVector3D<T> operator* (const BasicVector3D<T>& v, U a)
+{ return BasicVector3D<T>(v.x()*a, v.y()*a, v.z()*a); }
 
+//! Multiplication scalar by vector.
 //! @relates BasicVector3D
-//! We do not provide the operator form a*b:
-//! Though nice to write, and in some cases perfectly justified,
-//! in general it tends to make expressions more difficult to read.
-template <class T>
-inline T
-dotProduct(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
-{ return a.dot(b); }
+template <class T, class U>
+inline BasicVector3D<T> operator* (U a, const BasicVector3D<T>& v)
+{ return BasicVector3D<T>(a*v.x(), a*v.y(), a*v.z()); }
 
-//! Cross product.
+// vector*vector not supported
+//    (We do not provide the operator form a*b of the dot product:
+//     Though nice to write, and in some cases perfectly justified,
+//     in general it tends to make expressions more difficult to read.)
+
+//! Division vector by scalar.
 //! @relates BasicVector3D
-template <class T>
-inline BasicVector3D<T>
-crossProduct(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
-{ return a.cross(b); }
+template <class T, class U>
+inline BasicVector3D<T> operator/ (const BasicVector3D<T>& v, U a)
+{ return BasicVector3D<T>(v.x()/a, v.y()/a, v.z()/a); }
 
 //! Comparison of two vectors for equality.
 //! @relates BasicVector3D
 template <class T>
-inline bool
-operator==(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
+inline bool operator==(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
 { return (a.x()==b.x()&& a.y()==b.y()&& a.z()==b.z()); }
 
 //! Comparison of two vectors for inequality.
 //! @relates BasicVector3D
 template <class T>
-inline bool
-operator!=(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
+inline bool operator!=(const BasicVector3D<T>& a, const BasicVector3D<T>& b)
 { return (a.x()!=b.x() || a.y()!=b.y() || a.z()!=b.z()); }
 
+// -----------------------------------------------------------------------------
+// Quasi constructor
+// -----------------------------------------------------------------------------
+
+//! Creates a vector<double> as a wavevector with given wavelength and angles.
+//! Specifically needed for grazing-incidence scattering.
+BasicVector3D<double> vecOfLambdaAlphaPhi(const double _lambda, const double _alpha, const double _phi);
+
+
+// =============================================================================
+// ?? for API generation ??
+// =============================================================================
+
 template<> BA_CORE_API_ double BasicVector3D<double>::mag2() const;
 
 template<> BA_CORE_API_ double BasicVector3D<double>::mag() const;
@@ -317,6 +291,7 @@ template<> BA_CORE_API_ double BasicVector3D<double>::phi() const;
 
 template<> BA_CORE_API_ double BasicVector3D<double>::theta() const;
 
+//! \todo Replace by member function complex()
 BA_CORE_API_ BasicVector3D<std::complex<double> > toComplexVector(
         const BasicVector3D<double>& real_vector);
 
diff --git a/Core/Geometry/src/BasicVector3D.cpp b/Core/Geometry/src/BasicVector3D.cpp
index 82cb15e7ee8..aed8381b840 100644
--- a/Core/Geometry/src/BasicVector3D.cpp
+++ b/Core/Geometry/src/BasicVector3D.cpp
@@ -115,6 +115,14 @@ double BasicVector3D<double>::sin2Theta() const
     return mag2() == 0 ? 0 : magxy2()/mag2();
 }
 
+//! Returns this, trivially converted to complex type.
+template<>
+BasicVector3D<std::complex<double>> BasicVector3D<double>::complex() const
+{
+    return BasicVector3D<std::complex<double>>( v_[0], v_[1], v_[2] );
+}
+
+
 // -----------------------------------------------------------------------------
 // Combine two vectors
 // -----------------------------------------------------------------------------
@@ -171,4 +179,17 @@ BasicVector3D<std::complex<double> > toComplexVector(const BasicVector3D<double>
                                                 complex_t(real_vector.z()));
 }
 
+// -----------------------------------------------------------------------------
+// Quasi constructor
+// -----------------------------------------------------------------------------
+
+BasicVector3D<double> vecOfLambdaAlphaPhi(const double _lambda, const double _alpha, const double _phi)
+{
+    double k = PI2/_lambda;
+    return BasicVector3D<double>(
+        k*std::cos(_alpha) * std::cos(_phi),
+        -k*std::cos(_alpha) * std::sin(_phi),
+        k*std::sin(_alpha) );
+}
+
 }  // namespace Geometry
diff --git a/Core/Samples/src/Lattice.cpp b/Core/Samples/src/Lattice.cpp
index 9e6bc440dbf..d0abaddf333 100644
--- a/Core/Samples/src/Lattice.cpp
+++ b/Core/Samples/src/Lattice.cpp
@@ -75,7 +75,7 @@ void Lattice::initialize() const
 
 double Lattice::getVolume() const
 {
-    return std::abs(dotProduct(m_a1, crossProduct(m_a2, m_a3)));
+    return std::abs(m_a1.dot( m_a2.cross(m_a3)));
 }
 
 void Lattice::getReciprocalLatticeBasis(kvector_t& b1, kvector_t& b2,
@@ -144,12 +144,12 @@ Lattice Lattice::createTrigonalLattice(double a, double c)
 
 void Lattice::computeReciprocalVectors() const
 {
-    kvector_t a23 = crossProduct(m_a2, m_a3);
-    kvector_t a31 = crossProduct(m_a3, m_a1);
-    kvector_t a12 = crossProduct(m_a1, m_a2);
-    m_b1 = Units::PI2/dotProduct(m_a1, a23)*a23;
-    m_b2 = Units::PI2/dotProduct(m_a2, a31)*a31;
-    m_b3 = Units::PI2/dotProduct(m_a3, a12)*a12;
+    kvector_t a23 = m_a2.cross(m_a3);
+    kvector_t a31 = m_a3.cross(m_a1);
+    kvector_t a12 = m_a1.cross(m_a2);
+    m_b1 = Units::PI2/m_a1.dot(a23)*a23;
+    m_b2 = Units::PI2/m_a2.dot(a31)*a31;
+    m_b3 = Units::PI2/m_a3.dot(a12)*a12;
 }
 
 void Lattice::computeVectorsWithinRadius(const kvector_t& input_vector,
diff --git a/Core/Tools/inc/Bin.h b/Core/Tools/inc/Bin.h
index a8815bc9d9d..59c7e720cca 100644
--- a/Core/Tools/inc/Bin.h
+++ b/Core/Tools/inc/Bin.h
@@ -46,8 +46,8 @@ struct BA_CORE_API_ Bin1DKVector
     Bin1DKVector() : m_q_lower(), m_q_upper() {}
     Bin1DKVector(const kvector_t& lower, const kvector_t& upper)
         : m_q_lower(lower), m_q_upper(upper) {}
-    Bin1DKVector(double wavelength, const Bin1D& alpha_bin,
-                 const Bin1D& phi_bin);
+    Bin1DKVector(double wavelength, const Bin1D& alpha_bin, const Bin1D& phi_bin);
+
     kvector_t getMidPoint() const { return (m_q_lower + m_q_upper)/2.0; }
     kvector_t getDelta() const { return m_q_upper - m_q_lower; }
     kvector_t m_q_lower;  //!< lower bound of the bin
@@ -58,8 +58,8 @@ struct BA_CORE_API_ Bin1DKVector
 inline Bin1DKVector::Bin1DKVector(double wavelength, const Bin1D &alpha_bin, const Bin1D &phi_bin)
     : m_q_lower(), m_q_upper()
 {
-    m_q_lower.setLambdaAlphaPhi(wavelength, alpha_bin.m_lower, phi_bin.m_lower);
-    m_q_upper.setLambdaAlphaPhi(wavelength, alpha_bin.m_upper, phi_bin.m_upper);
+    m_q_lower = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_bin.m_lower, phi_bin.m_lower);
+    m_q_upper = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_bin.m_upper, phi_bin.m_upper);
 }
 
 //! @class Bin1DCVector
@@ -82,8 +82,8 @@ struct BA_CORE_API_ Bin1DCVector
 inline Bin1DCVector::Bin1DCVector(double wavelength, const Bin1D& alpha_bin, const Bin1D& phi_bin)
     : m_q_lower(), m_q_upper()
 {
-    m_q_lower.setLambdaAlphaPhi(wavelength, alpha_bin.m_lower, phi_bin.m_lower);
-    m_q_upper.setLambdaAlphaPhi(wavelength, alpha_bin.m_upper, phi_bin.m_upper);
+    m_q_lower = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_bin.m_lower, phi_bin.m_lower).complex();
+    m_q_upper = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_bin.m_upper, phi_bin.m_upper).complex();
 }
 
 #endif /* BIN_H_ */
-- 
GitLab