Sorry, I was wrong. I missed that z is in range -1..1. Original report is invalid, random.vonmisesvariate() always returns a value on the full circle.

However there is some inconsistency. For small kappa (<= 1e-6) result range is 0 to 2pi, for other kappa it is (mu%2pi)-pi to (mu%2pi)+pi. For consistency we should either shift a range for small kappa:

        if kappa <= 1e-6:
            return (mu % TWOPI) - _pi + TWOPI * random()

or normalize a result in another case:

        if u3 > 0.5:
            theta = (mu + _acos(f)) % TWOPI
        else:
            theta = (mu - _acos(f)) % TWOPI