How to Handle Two-Center Bipolar Coordinates?

In my problem, I want to integrate a $2D$ function $f(x,y)$ which explicitly depends on the vector $ \vec{r}_1=\vec{r}-\vec{R}_1 $ and $\vec{r}_2=\vec{r}-\vec{R}_2$, where $\vec{R}_1=(a,0)$ and $\vec{R}_2=(-a,0)$ are two fixed points. I found that the two-center bipolar coordinate maybe helpful which you can find some information about it here in Wikipedia.

However, I am not clear how to apply the integration exactly. The question is general (and maybe stupid): Since the two coordinate $r_1$ and $r_2$ are not orthogonal, would there be any problem when doing the change of variables $(x,y)$ to $(r_1,r_2)$?

Furthermore, how to express differential operators, like gradient, curl, in the new coordinates?

Any useful reference are welcome.

Solutions Collecting From Web of "How to Handle Two-Center Bipolar Coordinates?"