From 2a9db058c21a833e092d94555130f82c532f03ef Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Mon, 18 Jun 2018 12:13:03 +0200
Subject: [PATCH] Copy evaluation in IComputationTerm to function call operator

---
 .../GISASSpecularComputationTerm.cpp          | 18 +++++++++
 .../GISASSpecularComputationTerm.h            |  2 +
 Core/Computation/IComputationTerm.h           |  3 ++
 .../Computation/ParticleLayoutComputation.cpp | 14 +++++++
 Core/Computation/ParticleLayoutComputation.h  |  2 +
 .../RoughMultiLayerComputation.cpp            | 39 +++++++++++++++++++
 Core/Computation/RoughMultiLayerComputation.h |  2 +
 7 files changed, 80 insertions(+)

diff --git a/Core/Computation/GISASSpecularComputationTerm.cpp b/Core/Computation/GISASSpecularComputationTerm.cpp
index e0ad9a7fb4e..aff75091f5f 100644
--- a/Core/Computation/GISASSpecularComputationTerm.cpp
+++ b/Core/Computation/GISASSpecularComputationTerm.cpp
@@ -35,6 +35,24 @@ void GISASSpecularComputationTerm::eval(
             evalSingle(it);
 }
 
+void GISASSpecularComputationTerm::operator()(SimulationElement& elem) const
+{
+    if (mp_multilayer->requiresMatrixRTCoefficients())
+        return;
+
+    if (elem.isSpecular()) {
+        complex_t R = mp_fresnel_map->getInCoefficients(elem, 0)->getScalarR();
+        double sin_alpha_i = std::abs(std::sin(elem.getAlphaI()));
+        if (sin_alpha_i == 0.0)
+            sin_alpha_i = 1.0;
+        const double solid_angle = elem.getSolidAngle();
+        if (solid_angle <= 0.0)
+            return;
+        const double intensity = std::norm(R) * sin_alpha_i / solid_angle;
+        elem.setIntensity(intensity);
+    }
+}
+
 void GISASSpecularComputationTerm::evalSingle(const std::vector<SimulationElement>::iterator& iter) const
 {
     complex_t R = mp_fresnel_map->getInCoefficients(*iter, 0)->getScalarR();
diff --git a/Core/Computation/GISASSpecularComputationTerm.h b/Core/Computation/GISASSpecularComputationTerm.h
index 2c385ec3ea5..ab0e38507c1 100644
--- a/Core/Computation/GISASSpecularComputationTerm.h
+++ b/Core/Computation/GISASSpecularComputationTerm.h
@@ -30,6 +30,8 @@ public:
     void eval(ProgressHandler* progress, const std::vector<SimulationElement>::iterator& begin_it,
               const std::vector<SimulationElement>::iterator& end_it) const override;
 
+    void operator()(SimulationElement& elem) const override;
+
 private:
     void evalSingle(const std::vector<SimulationElement>::iterator& iter) const;
 };
diff --git a/Core/Computation/IComputationTerm.h b/Core/Computation/IComputationTerm.h
index cb7479d9243..25ee23060e9 100644
--- a/Core/Computation/IComputationTerm.h
+++ b/Core/Computation/IComputationTerm.h
@@ -44,6 +44,9 @@ public:
                       const std::vector<SimulationElement>::iterator& begin_it,
                       const std::vector<SimulationElement>::iterator& end_it) const =0;
 
+    //! Calculate scattering intensity for the given SimulationElement
+    virtual void operator()(SimulationElement& elem) const=0;
+
     //! Merges its region map into the given one (notice non-const reference parameter)
     void mergeRegionMap(std::map<size_t, std::vector<HomogeneousRegion>>& region_map) const;
 
diff --git a/Core/Computation/ParticleLayoutComputation.cpp b/Core/Computation/ParticleLayoutComputation.cpp
index 2967678b8d9..6f4288cd590 100644
--- a/Core/Computation/ParticleLayoutComputation.cpp
+++ b/Core/Computation/ParticleLayoutComputation.cpp
@@ -60,3 +60,17 @@ void ParticleLayoutComputation::eval(ProgressHandler* progress,
         }
     }
 }
+
+void ParticleLayoutComputation::operator()(SimulationElement& elem) const
+{
+    double alpha_f = elem.getAlphaMean();
+    size_t n_layers = mp_multilayer->numberOfLayers();
+    if (n_layers > 1 && alpha_f < 0) {
+        return; // zero for transmission with multilayers (n>1)
+    } else {
+        elem.addIntensity(mP_strategy->evaluate(elem) * m_surface_density);
+    }
+    if (mP_progress_counter) {
+        mP_progress_counter->stepProgress();
+    }
+}
diff --git a/Core/Computation/ParticleLayoutComputation.h b/Core/Computation/ParticleLayoutComputation.h
index e33f326251e..7b9201c3e84 100644
--- a/Core/Computation/ParticleLayoutComputation.h
+++ b/Core/Computation/ParticleLayoutComputation.h
@@ -41,6 +41,8 @@ public:
               const std::vector<SimulationElement>::iterator& begin_it,
               const std::vector<SimulationElement>::iterator& end_it) const override;
 
+    void operator()(SimulationElement& elem) const override;
+
 private:
     std::unique_ptr<const IInterferenceFunctionStrategy> mP_strategy;
     double m_surface_density;
diff --git a/Core/Computation/RoughMultiLayerComputation.cpp b/Core/Computation/RoughMultiLayerComputation.cpp
index 53ed64d310b..03eb5bbb351 100644
--- a/Core/Computation/RoughMultiLayerComputation.cpp
+++ b/Core/Computation/RoughMultiLayerComputation.cpp
@@ -64,6 +64,45 @@ void RoughMultiLayerComputation::eval(ProgressHandler* progress,
     }
 }
 
+void RoughMultiLayerComputation::operator()(SimulationElement& elem) const
+{
+    if (elem.getAlphaMean()<0.0)
+        return;
+    kvector_t q = elem.getMeanQ();
+    double wavelength = elem.getWavelength();
+    double autocorr(0.0);
+    complex_t crosscorr(0.0, 0.0);
+
+    std::vector<complex_t > rterm( mp_multilayer->numberOfLayers()-1 );
+    std::vector<complex_t > sterm( mp_multilayer->numberOfLayers()-1 );
+
+    for (size_t i=0; i<mp_multilayer->numberOfLayers()-1; i++){
+        rterm[i] = get_refractive_term(i, wavelength);
+        sterm[i] = get_sum8terms(i, elem);
+    }
+    for (size_t i=0; i<mp_multilayer->numberOfLayers()-1; i++) {
+        const LayerRoughness *rough =
+            mp_multilayer->layerBottomInterface(i)->getRoughness();
+        if(rough)
+            autocorr += std::norm( rterm[i] ) * std::norm( sterm[i] ) * rough->getSpectralFun(q);
+    }
+    // cross correlation between layers
+    if (mp_multilayer->crossCorrLength() != 0.0) {
+        for(size_t j=0; j<mp_multilayer->numberOfLayers()-1; j++){
+            for(size_t k=0; k<mp_multilayer->numberOfLayers()-1; k++) {
+                if(j==k) continue;
+                crosscorr += rterm[j]*sterm[j]*
+                    mp_multilayer->crossCorrSpectralFun(q,j,k)*
+                    std::conj(rterm[k])*std::conj(sterm[k]);
+            }
+        }
+    }
+    //! @TODO clarify complex vs double
+    elem.addIntensity((autocorr+crosscorr.real())*M_PI/4./wavelength/wavelength);
+
+    mP_progress_counter->stepProgress();
+}
+
 double RoughMultiLayerComputation::evaluate(const SimulationElement& sim_element) const
 {
     if (sim_element.getAlphaMean()<0.0)
diff --git a/Core/Computation/RoughMultiLayerComputation.h b/Core/Computation/RoughMultiLayerComputation.h
index 2ea9914ece6..49162dbe17b 100644
--- a/Core/Computation/RoughMultiLayerComputation.h
+++ b/Core/Computation/RoughMultiLayerComputation.h
@@ -39,6 +39,8 @@ public:
               const std::vector<SimulationElement>::iterator& begin_it,
               const std::vector<SimulationElement>::iterator& end_it) const override;
 
+    void operator()(SimulationElement& elem) const override;
+
 private:
     double evaluate(const SimulationElement& sim_element) const;
     complex_t get_refractive_term(size_t ilayer, double wavelength) const;
-- 
GitLab