def draw_bundle(rng, flux):
    """
    Generate a bundle impact point, zenith angle, and multiplicity proportional
    to the total flux
    """
    max_flux = float(flux(1.5, 1., flux.min_multiplicity))
    while True:
        depth = rng.uniform(1.5, 2.5)
        ct = rng.uniform(0, 1)
        m = flux.min_multiplicity + rng.integer(flux.max_multiplicity-flux.min_multiplicity)
        if rng.uniform(0, max_flux) < float(flux(depth, ct, m)):
            return depth, ct, m