I kind of follow your argument, but the implementation in this case is above my pay grade.

My reply to the previous response supplies more information. The weighting function is always real and has no poles and is well behaved. However, Mathematica often tells me that the convergence is too slow or that the integrand may be oscillating....so I used one or another of the methods in Mathematica to speed things up. I also play with the accuracy and precision setting. At the same settings, the results as described above differ between the two coordinate systems.

The integrals, at first blush, resemble elliptical integrals but are more complex. There are no analytical solutions, so the calculations must be made numerically.

Well, I hate to say it, but it all depends.

What I have is a function that depends on the vector between volume elements within two volumes. I am integrating this function over both volumes. The volumes are spheres. I can define the function or the volumes in spherical or cylindrical coordinates. The function essentially "weights" the vector for each volume element pair.

When I set the function to unity, the two integrals produce the same volumes. When I allow the function to take on a value related to the length of the vector between volume elements, the value of this integral differs as the ratio of the sphere radii is varied. Hence, I get a range of values over the range of the radii ratio. When I switch from cylindrical to spherical coordinates, the value of the integral is not the same between the two coordinate systems for a given ratio. Hence, I get a type of scatter plot, but the scatter plot is different between the two coordinate systems.