Understanding decomposition tracking

decomposition_tracking.cpp

void decomposition_tracking(const Ray& ray, const Material& material, float& t_out, float3& transmittance) { 
    auto const cm = material.control_medium();

    float const rc  = rng_.random_float();
    float const t_c = ray.min_t - std::log(1.f - rc) / cm.minorant_mu_t;

    if (t_c > ray.max_t) {
        t_out = ray.max_t;
        transmittance = float3(1.f);
        return;
    }
    
    float const mt = cm.majorant_mu_t - cm.minorant_mu_t;

    auto mu = material.collision_coefficients();
    mu.a -= cm.minorant_mu_a;
    mu.s -= cm.minorant_mu_s;

    float const mu_n = cm.majorant_mu_t - cm.minorant_mu_t - mt;

    for (float t_r = ray.min_t;;) {
        float const r0 = rng_.random_float();
        t_r -= std::log(1.f - r0) / mt;

        if (t > t_c) {
            float const r1 = rng_.random_float();

            if (r1 < cm.minorant_mu_a / cm.minorant_mu_t) {
                t_out = ray.max_t;
                transmittance = float3(0.f);
                return;
            } else {
                t_out = t_c;
                transmittance = float3(1.f);
                return;
            }
        } else {
            float const r1 = rng_.random_float();

            if (r1 < mu.a / mt) {
                t_out = ray.max_t;
                transmittance = float3(0.f);
                return;
            } else if (r1 < 1.f - mu_n / mt) {
                t_out = t_r;
                transmittance = float3(1.f);
                return;
            }
        }
    } 
}