Support source and sink below horizon for scalar scattering

Resample/Specular/ComputeFluxMagnetic.cpp

@@ -174,6 +174,9 @@ std::vector<MatrixFlux> computeFlux(const SliceStack& slices, const std::vector<
 
 Fluxes Compute::SpecularMagnetic::fluxes(const SliceStack& slices, const kvector_t& k, bool forward)
 {
+    if (slices.size()>1 && k.z()>0)
+        throw std::runtime_error(
+            "source or detector below horizon not yet implemented for polarized scattering");
     std::vector<complex_t> kz = Compute::Kz::computeReducedKz(slices, k);
     ASSERT(slices.size() == kz.size());
 

Resample/Specular/ComputeFluxScalar.cpp

@@ -68,7 +68,7 @@ std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double
 }
 
 std::vector<Eigen::Vector2cd> computeTR(const SliceStack& slices, const std::vector<complex_t>& kz,
-                                        const RoughnessModel& r_model)
+                                        const RoughnessModel& r_model, bool top_exit)
 {
     const size_t N = slices.size();
     std::vector<Eigen::Vector2cd> TR(N, {1., 0.});
@@ -77,37 +77,45 @@ std::vector<Eigen::Vector2cd> computeTR(const SliceStack& slices, const std::vec
         return TR;
 
     if (kz[0] == 0.0) { // If kz in layer 0 is zero, R0 = -T0 and all others equal to 0
+        ASSERT(top_exit);
         TR[0] = {1.0, -1.0};
         for (size_t i = 1; i < N; ++i)
             TR[i].setZero();
         return TR;
     }
+    ASSERT(! (!top_exit && kz[N-1] == 0.0) ); // forbidden as case top_exit=true includes k_z=0
+
+    // Index reversal for bottom exit
+    std::vector<size_t> X(N, 0);
+    for (size_t i=0; i<N; ++i)
+        X[i] = top_exit ? i : N-1-i;
 
     // Calculate transmission/refraction coefficients t_r for each layer, from bottom to top.
-    TR.back() = {1.0, 0.0};
+    TR[X[N-1]] = {1.0, 0.0};
     std::vector<complex_t> factors(N - 1);
     for (int i = N - 2; i >= 0; i--) {
-        double sigma = 0.0;
-        if (const auto roughness = slices.bottomRoughness(i))
-            sigma = roughness->getSigma();
+        size_t jthis = X[i];
+        size_t jlast = X[i + 1];
+        const auto roughness = slices.bottomRoughness(jthis); // TODO verify
+        const double sigma = roughness ? roughness->getSigma() : 0.;
 
-        const auto [mp, mm] = transition(kz[i], kz[i + 1], sigma, r_model);
+        const auto [mp, mm] = transition(kz[jthis], kz[jlast], sigma, r_model);
 
-        const complex_t delta = exp_I(kz[i] * slices[i].thicknessOr0());
+        const complex_t delta = exp_I(kz[jthis] * slices[jthis].thicknessOr0());
 
-        complex_t S = delta / (mp + mm * TR[i + 1](1));
+        complex_t S = delta / (mp + mm * TR[jlast](1));
         factors[i] = S;
-        TR[i](1) = delta * (mm + mp * TR[i + 1](1)) * S;
+        TR[jthis](1) = delta * (mm + mp * TR[jlast](1)) * S;
     }
 
     // Now correct all amplitudes by dividing by the remaining factors in forward direction.
     // At some point this divison underflows; from there on all further amplitudes are set to zero.
     complex_t dampingFactor = 1;
-    for (size_t j = 1; j < N; ++j) {
-        dampingFactor = dampingFactor * factors[j - 1];
+    for (size_t i = 1; i < N; ++i) {
+        dampingFactor = dampingFactor * factors[i - 1];
 
-        TR[j](0) = dampingFactor;
-        TR[j](1) *= dampingFactor;
+        TR[X[i]](0) = dampingFactor;
+        TR[X[i]](1) *= dampingFactor;
     }
 
     return TR;
@@ -121,10 +129,12 @@ std::vector<Eigen::Vector2cd> computeTR(const SliceStack& slices, const std::vec
 
 Fluxes Compute::SpecularScalar::fluxes(const SliceStack& slices, const kvector_t& k)
 {
+    const bool top_exit = k.z() <= 0; // true if source or detector pixel are at z>=0
     std::vector<complex_t> kz = Compute::Kz::computeReducedKz(slices, k);
     ASSERT(slices.size() == kz.size());
 
-    const std::vector<Eigen::Vector2cd> TR = computeTR(slices, kz, slices.roughnessModel());
+    const std::vector<Eigen::Vector2cd> TR = computeTR(
+        slices, kz, slices.roughnessModel(), top_exit);
 
     Fluxes result;
     for (size_t i = 0; i < kz.size(); ++i)