Integrating $\mathbf{z}$ over a function of shifted vectors $\|\mathbf{z} - \mathbf{x}\|$ and $\|\mathbf{z} - \mathbf{y}\|$ Fix two vectors $\mathbf{x} = \left[ \begin{matrix} x_{1} \\ x_{2} \end{matrix} \right]$ and $\mathbf{y} = \left[ \begin{matrix} y_{1} \\ y_{2} \end{matrix} \right]$. Let's suppose I am interested in the following integral:
$$
\int_{-\infty}^{\infty} dz_{1} \int_{-\infty}^{\infty} dz_{2} \ f\left( \| \mathbf{z} - \mathbf{x} \|^{2}, \| \mathbf{z} - \mathbf{y} \|^{2} \right)
$$
Where I am integrating over all $\mathbf{z} = \left[ \begin{matrix} z_{1} \\ z_{2} \end{matrix} \right] \in \mathbb{R}^{2}$, and $f$ is a function of the following:
$$
\| \mathbf{z} - \mathbf{x} \|^{2} = ( z_{1} - x_{1} )^{2} + ( z_{2} - x_{2} )^{2} \\
\| \mathbf{z} - \mathbf{y} \|^{2} = ( z_{1} - y_{1} )^{2} + ( z_{2} - y_{2} )^{2}
$$
Is there a change of coordinates that would make this computation a little simpler? Something a long the lines of polar coordinates would be nice.
I have tried defining the variables $U = \| \mathbf{z} - \mathbf{x} \|^{2}$ and $V = \| \mathbf{z} - \mathbf{y} \|^{2}$, but I don't think that these variables work well because I don't know how to define the inverse transformation $z_{1} = z_{1}(U,V)$ and $z_{2} = z_{2}(U,V)$ (from which to get the Jacobian matrix).
I have also thought about using elliptical coordinates enter link description here, but can't seem to understand how to fit this to use in my problem (should I rotate and shift my vectors $\mathbf{x}$ and $\mathbf{y}$ so that I match the definition used in the wikipedia link, in the sense that $a = \|\mathbf{x} - \mathbf{y}\|$?).
Thanks in advance for any help.
 A: You might be interested in the two center bipolar coordinate system (https://en.wikipedia.org/wiki/Two-center_bipolar_coordinates).  Here the independent variables are exactly the relevant distances, if you choose x and y as the points through which the line defining this coordinate system pass.
A: So if I define the following two variables:
$$
U := \sqrt{ ( z_{1} - x_{1} )^2 + ( z_{2} - x_{2} )^2 } \\
V := \sqrt{ ( z_{1} - y_{1} )^2 + ( z_{2} - y_{2} )^2 }
$$
Then I think my new variables should run over the limits of integration $U,V \in (0,\infty)$. Also I find that the Jacobian determinant is given by:
$$
\left| \frac{\partial \left( z_{1}, z_{2} \right)}{ \partial \left( U, V \right) } \right| \ = \ \frac{ 2 U V }{ \sqrt{ - \left( (U+V)^2-d^2 \right) \left( (U-V)^2-d^2 \right) } }
$$
...where I have defined the constant $d = ||\mathbf{x} - \mathbf{y}|| = \sqrt{( x_{1} - y_{1} )^2 + ( x_{2} - y_{2} )^2}$.
Therefore, my integral of interest now becomes written as:
$$
\int_{-\infty}^{\infty} dz_{1} \int_{-\infty}^{\infty} dz_{2} \ f\left( || \mathbf{z} - \mathbf{x} ||^2, || \mathbf{z} - \mathbf{y} ||^2 \right) \ = \ \int_{0}^{\infty} dU \int_{0}^{\infty} dV \ f\left( U^2, V^2 \right) \frac{ 2 U V }{ \sqrt{ - \left( (U+V)^2-d^2 \right) \left( (U-V)^2-d^2 \right) } }
$$ 
My issue now is that the Jacobian is complex-valued for large enough $X$ and $Y$! What can this possibly mean? Is this telling me that these are not valid coordinates?
