store ratio's in separate array

Resample/Specular/ComputeFluxScalar.cpp

@@ -85,7 +85,8 @@ std::vector<Spinor> computeTR(const SliceStack& slices, const std::vector<comple
     }
 
     // Parratt algorithm, step 1: compute r/t ratios from bottom to top.
-    std::vector<complex_t> tfactor(N - 1);
+    std::vector<complex_t> tfactor(N - 1); // transmission damping t_{i+1} / t_{i}
+    std::vector<complex_t> ratio(N, 0.);   // Parratt's x = r/t
     for (size_t i = N - 1; i > 0; i--) {
         const size_t jthis = X[i - 1];
         const size_t jlast = X[i];
@@ -95,16 +96,17 @@ std::vector<Spinor> computeTR(const SliceStack& slices, const std::vector<comple
         const auto [slp, slm] = transition(kz[jthis], kz[jlast], sigma, r_model);
 
         const complex_t delta = exp_I(kz[jthis] * slices[jthis].thicknessOr0());
-        const complex_t f = delta / (slp + slm * TR[jlast].v);
+        const complex_t f = delta / (slp + slm * ratio[jlast]);
         tfactor[i - 1] = f;
-        TR[jthis].v = delta * (slm + slp * TR[jlast].v) * f;
+        ratio[jthis] = delta * (slm + slp * ratio[jlast]) * f;
+        TR[jthis].v = ratio[jthis];
     }
 
     // Parrat algorithm, step 2: compute r and t from top to bottom.
     complex_t t = 1;
     for (size_t i = 1; i < N; ++i) {
         t = t * tfactor[i - 1];
-        TR[X[i]] = Spinor(t, TR[X[i]].v * t);
+        TR[X[i]] = Spinor(t, ratio[X[i]] * t);
     }
 
     return TR;