optimized replication through mediation

Sim/Computation/RoughMultiLayerContribution.cpp

@@ -100,34 +100,28 @@ complex_t get_sum8terms(const ReSample& re_sample, size_t i_layer, const Diffuse
     return term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8;
 }
 
-// Fourier transform of the correlation function of roughnesses between the interfaces
-double crossCorrSpectralFun(const R3& qvec, const SliceStack& stack, size_t j, size_t k)
+double underlyingInterfaceReplication(const R3& qvec, const SliceStack& stack, size_t j)
 {
-    ASSERT(j < k);
-
-    const double distance = std::abs(stack[j].hig() - stack[k].hig());
-    const AutocorrelationModel* rough_j = stack[j].topRoughness()->autocorrelationModel();
-    const AutocorrelationModel* rough_k = stack[k].topRoughness()->autocorrelationModel();
-    const double sigma_j = rough_j->sigma();
-    const double sigma_k = rough_k->sigma();
+    const CrosscorrelationModel* crosscorr_j = stack[j].topRoughness()->crosscorrelationModel();
 
-    if (sigma_j <= 0 || sigma_k <= 0)
+    if (!crosscorr_j)
         return 0.0;
 
-    const CrosscorrelationModel* crosscorr_j = stack[j].topRoughness()->crosscorrelationModel();
+    const double distance = std::abs(stack[j].hig() - stack[j + 1].hig());
+    return crosscorr_j->replicationFactor(qvec, distance);
+}
 
-    if (!crosscorr_j)
+double crossInterfaceSpectralFun(double spectrum_j, double spectrum_k, double sigma_j,
+                                 double sigma_k)
+{
+    if (sigma_j <= 0 || sigma_k <= 0)
         return 0.0;
 
-    return 0.5
-           * ((sigma_k / sigma_j) * rough_j->spectralFunction(qvec)
-              + (sigma_j / sigma_k) * rough_k->spectralFunction(qvec))
-           * crosscorr_j->replicationFactor(qvec, distance);
+    return 0.5 * ((sigma_k / sigma_j) * spectrum_j + (sigma_j / sigma_k) * spectrum_k);
 }
 
 } // namespace
 
-
 double Compute::roughMultiLayerContribution(const ReSample& re_sample, const DiffuseElement& ele)
 {
     if (ele.alphaMean() < 0.0)
@@ -155,19 +149,35 @@ double Compute::roughMultiLayerContribution(const ReSample& re_sample, const Dif
         }
 
     const size_t n_interfaces = roughStack.size();
+    std::vector<double> spectrum(n_interfaces, 0);
 
     // auto correlation in each layer (first term in final expression in Eq (3) of Schlomka et al)
-    for (size_t i = 0; i < n_interfaces; i++)
-        autocorr += std::norm(rterm[i] * sterm[i])
-                    * roughStack[i].topRoughness()->autocorrelationModel()->spectralFunction(q);
+    for (size_t i = 0; i < n_interfaces; i++) {
+        spectrum[i] = roughStack[i].topRoughness()->autocorrelationModel()->spectralFunction(q);
+        autocorr += std::norm(rterm[i] * sterm[i]) * spectrum[i];
+    }
 
     // cross correlation between layers (second term in loc. cit.)
-    for (size_t j = 0; j < n_interfaces; j++) {
-        if (!roughStack[j].topRoughness()->crosscorrelationModel())
+    std::vector<double> underlying_interface_replication(n_interfaces, 0); // the last is unused
+    std::vector<double> sigma(n_interfaces);
+
+    for (size_t j = 0; j < n_interfaces - 1; j++) {
+        sigma[j] = roughStack[j].topRoughness()->sigma();
+        underlying_interface_replication[j] = underlyingInterfaceReplication(q, roughStack, j);
+    }
+
+    for (size_t j = 0; j < n_interfaces - 1; j++) {
+        if (underlying_interface_replication[j] == 0)
             continue;
-        for (size_t k = j + 1; k < n_interfaces; k++)
+        double accumulated_replication = underlying_interface_replication[j];
+        for (size_t k = j + 1; k < n_interfaces; k++) {
             crosscorr += (rterm[j] * sterm[j] * std::conj(rterm[k] * sterm[k])).real()
-                         * ::crossCorrSpectralFun(q, roughStack, j, k);
+                         * crossInterfaceSpectralFun(spectrum[j], spectrum[k], sigma[j], sigma[k])
+                         * accumulated_replication;
+            if (underlying_interface_replication[k] == 0)
+                break;
+            accumulated_replication *= underlying_interface_replication[k];
+        }
     }
 
     const double k0 = (2 * pi) / wavelength;