This decomposition relation motivates to search for an algorithm that is capable of faithfully constructing the individual processes represented by the rhs of Eq. (REF ), upon sampling from a free-space density functions only.
In case of a reflective boundary, sampling according to the renormalized free-space density corresponds to the first two factors on the rhs. We now detail an iterative bisection method to faithfully model the FP time process described by the third factor on the rhs. A similar construction has been described in [1]} to determine the last reflection time for a 1D problem in a different context, but it works equally well in 2 and 3D for our problem at hand.
