Copy evaluation in IComputationTerm to function call operator

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();

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;
 };

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;
 

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();
+    }
+}

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;

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)

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;