I believe my formula is right, also expression. |
Both integral limits are wrong.
The outer integral has to run from 0 to 80 (R+r) and the the upper limit of the inner integral gets imaginary for values greater than 40 (which is the reason for the error). So replace it by its real part which sets it to zero for values greater than 40 and this happens to be exactyl what is needed. Furthermore you should either use the second value of z or switch the integral limits, but thats just a matter of sign.
You formula will only calculate the volume of the upper half of a quarter of the donut, so you have to multiply it by 8.
If you let the outer integral go from -80 to 80 you get the upper half (you have chosen the positive root for z) of the right half (y>0) torus and have to multiply by 4 therefore.
For better accuracy you should decrease the value of TOL - default is 10^-3 which gets you a slighty different result compared to the "exact" v.boia.

There is no symbolic evaluation necessary:

And if you like to use units, too, you have to redifine z(x,y) - best done in a general form using R and r as arguments:
