Integrating for the volume of a torus in SciPy

[Nitram]

My guess is that for z_low and z_high it should be lambda theta, rho:... because the documentation says qfun(x, y) and rfun(x, y)

It turns out you are correct! Thanks!
I ran the code

import numpy as np
from scipy.integrate import tplquad

R = 4
r = 1

# Volume integral
# Order of integration is z, rho, theta
f1 = lambda z, rho, theta: 2 * rho
# set theta limits
a, b = 0, 2 * np.pi
# set rho limits
rho_low = lambda rho: R - r
rho_high = lambda rho: R + r
# set z limits
z_low = lambda theta, rho: 0
z_high = lambda theta, rho: np.sqrt(r**2 - (rho - R)**2)
# calculate integral
V, _ = tplquad(f1, a, b, rho_low, rho_high, z_low, z_high)
print(V)

which gives the correct result. I did not have to integrate in the first octant only and then multiply by 8 after all . How come the order 'lambda rho, theta' did not work and 'lambda theta, rho' does?