How can I plot the gradient field in bipolar coordinates?

Prelude: bipolar coordinates

I use $(u, v)$ rather than $(\sigma, \tau)$ throughout. The transformation is given by$$\begin{align*} x &= \frac{\sinh v}{\cosh v - \cos u} \\ y &= \frac{\sin u}{\cosh v - \cos u}.\end{align*}$$

Not so well known is the inverse transformation (which we will need for plotting):$$\begin{align*} u &= \tan^{-1} \frac{2y}{x^2 + y^2 - 1} \\ v &= \tanh^{-1} \frac{2x}{x^2 + y^2 + 1}.\end{align*}$$

The scale factors (Lamé coefficients) for both coordinates are the same:$$ h_u = h_v = h = \frac{1}{\cosh v - \cos u}.$$

The local orthonormal basis is related to the standard Cartesian basis according to$$\begin{align*} \mathbf{a}_u &= h (-S \,\mathbf{a}_x + C \,\mathbf{a}_y) \\ \mathbf{a}_v &= h (-C \,\mathbf{a}_x - S \,\mathbf{a}_y), \\\end{align*}$$ where $$\begin{align*} C &= \cos u \cosh v - 1 \\ S &= \sin u \sinh v.\end{align*}$$

Implementing this in Mathematica:

(* Coordinate transformations *) xBipolar[u_, v_] := Sinh[v] / (Cosh[v] - Cos[u]); yBipolar[u_, v_] := Sin[u] / (Cosh[v] - Cos[u]); (* Inverse coordinate transformations *) uBipolar[x_, y_] := ArcTan[x^2 + y^2 - 1, 2 y]; vBipolar[x_, y_] := ArcTanh[2 x / (x^2 + y^2 + 1)]; uvBipolar[x_, y_] := {uBipolar, vBipolar} @@ {x, y} // Through // Evaluate; (* Scale factors (both are the same) *) hBipolar[u_, v_] := 1 / (Cosh[v] - Cos[u]); (* Abbreviations *) cBipolar[u_, v_] := Cos[u] Cosh[v] - 1; sBipolar[u_, v_] := Sin[u] Sinh[v]; (* Cartesian components of local orthonormal basis *) uVectorBipolar[u_, v_] := hBipolar[u, v] {-sBipolar[u, v], cBipolar[u, v]} // Evaluate; vVectorBipolar[u_, v_] := hBipolar[u, v] {-cBipolar[u, v], -sBipolar[u, v]} // Evaluate; 

Function

I have taken $v_0 = 1$. We have$$ \frac{V}{V_0} = \sum_{n = 0}^\infty \frac{1}{\pi} \sqrt{2 (\cosh v - \cos u)} \cdot \lambda_n \cdot \frac{Q_{n-1/2}(\cosh v_0)}{P_{n-1/2}(\cosh v_0)} \cdot P_{n-1/2}(\cosh v) \cos (n u).$$

(* lambda-bar *) lambdaBar[0] = 1; lambdaBar[n_] /; n > 0 = 2; (* Expansion terms *) v0 = 1; term[n_][u_, v_] := ( 1 / Pi Sqrt[2 (Cosh[v] - Cos[u])] lambdaBar[n] LegendreQ[n - 1/2, Cosh[v0]] / LegendreP[n - 1/2, Cosh[v0]] LegendreP[n - 1/2, Cosh[v]] Cos[n u] ); (* Partial sum for V/V_0 *) partialSum[nMax_][u_, v_] := Sum[term[n][u, v], {n, 0, nMax}]; 

As mentioned in the comments, the LegendreQ factor isn't real, so we take the real part.

First we figure out how many terms we need to plot; I am guessing that the boundary condition which gives rise to the Fourier series is $V/V_0 = 1$ along the circle $v = v_0$:

nMaxValues = {0, 1, 2, 5}; Plot[ Table[ partialSum[nMax][u, v0] // Re , {nMax, nMaxValues} ] // Evaluate , {u, 0, 2 Pi} , PlotLegends -> LineLegend[nMaxValues, LegendLabel -> "nMax"] ] 

We see that 5 terms are enough. We certainly don't need 150.

(* Assume function intended for v < v_0 only *) regionFun = Function[{x, y}, Abs[vBipolar[x, y]] < v0]; (* Plot V/V_0 *) Plot3D[ partialSum[5] @@ uvBipolar[x, y] // Re // Evaluate , {x, -3, 3}, {y, -3, 3} , Exclusions -> None , RegionFunction -> regionFun ] 

Gradient

Next we take the gradient. Actually taking the derivative of a Fourier series is non-trivial. The $u$-derivative of $\cos(nu)$ introduces an extra factor of $n$, and if the coefficients do not go to zero fast enough, the term-by-term derivative will not converge. In this case though, the coefficients do go to zero fast enough for us to take a term-by-term derivative.

First define $\partial V / {\partial u}$ and $\partial V / {\partial v}$:

(* Derivatives of expansion terms *) termUDerivative[n_][u_, v_] := D[term[n][u, v], u] // Evaluate; termVDerivative[n_][u_, v_] := D[term[n][u, v], v] // Evaluate; (* Derivatives partial sum *) partialSumUDerivative[nMax_][u_, v_] := Sum[termUDerivative[n][u, v], {n, 0, nMax}]; partialSumVDerivative[nMax_][u_, v_] := Sum[termVDerivative[n][u, v], {n, 0, nMax}]; 

Since the scale factors for both coordinates are equal, the gradient is given by $$ \nabla V = \frac{1}{h} \left( \frac{\partial V}{\partial u} \,\mathbf{a}_u + \frac{\partial V}{\partial v} \,\mathbf{a}_v \right):$$

(* Partial sum for gradient of V/V_0 *) gradientPartialSum[nMax_][u_, v_] := 1 / hBipolar[u, v] * Plus[ partialSumUDerivative[nMax][u, v] uVectorBipolar[u, v], partialSumVDerivative[nMax][u, v] vVectorBipolar[u, v] ]; 

Finally we can plot the gradient (I chose StreamPlot because VectorPlot arrows are too small):

Show[ ContourPlot[ partialSum[5] @@ uvBipolar[x, y] // Re // Evaluate , {x, -3, 3}, {y, -3, 3} , AspectRatio -> Automatic , ContourShading -> None , Exclusions -> None , RegionFunction -> regionFun ], StreamPlot[ gradientPartialSum[5] @@ uvBipolar[x, y] // Re // Evaluate , {x, -3, 3}, {y, -3, 3} , RegionFunction -> regionFun ] ]