Plotting a vector field from spherical coordinates

Question

I am trying to plot an integral that is in spherical coordinates, but I am a bit lost. I think that my only issue is converting to cartesian, everything I have seen on how to do this has gone a bit over my head. Here is my code (sorry about the formatting):

$$ \vec{E}=\sigma k \int_0^R\int_0^{2\pi}[\frac{(r\sin\theta\cos\phi-r'\cos\phi')\textbf{i} +(r\sin\theta\sin\phi-r'\sin\phi')\textbf{j}+(r\cos\theta)\textbf{k}}{(r^2+r'^2-2rr'\cos\phi\cos\phi'\sin\theta-2rr'\sin\theta\sin\phi\sin\phi')^{3/2}}]r'd\phi'dr' $$

\[ScriptR][ϕ_, r_, θ_, ϕ1_, 
   r1_] := {r*Sin[θ]*Cos[ϕ] - r1*Cos[ϕ1], 
   r*Sin[θ]*Sin[ϕ] - r1*Sin[ϕ1], r*Cos[θ]};

\[ScriptR]Norm[ϕ_, r_, θ_, ϕ1_, r1_] := 
  Sqrt[(r*Sin[θ]*Cos[ϕ] - 
        r1*Cos[ϕ1])^2 + (r*Sin[θ]*Sin[ϕ] - 
        r1*Sin[ϕ1])^2 + (r*Cos[θ])^2] // Simplify;

ele[ϕ_?NumericQ, r_?NumericQ, θ_?NumericQ] := σ*k*
   NIntegrate[(\[ScriptR][ϕ, r, θ, ϕ1, 
        r1]/\[ScriptR]Norm[ϕ, r, θ, ϕ1, r1]^3)*
     r1, {ϕ1, 0, 2*π}, {r1, 0, R}];

VectorPlot3D[
 ele[ϕ, r, θ] /. {σ -> 200, k -> 9*10^9, R -> 0.5, 
   r*Sin[θ]*Cos[ϕ] -> x, r*Sin[θ]*Sin[ϕ] -> y,
    r*Cos[θ] -> z}, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]

I am worried part of my issue is in the function that I made, but I am most certain that calling replace all like this is not a viable way to transform coordinates.

Edit: added the integral I am trying to solve. The \[ScriptR] is the vector in the numerator, and \[ScriptR]Norm is the cube root of the denominator

Answer 1

VectorPlot3D expects a function of Cartesian $(x,y,z)$ coordinates, but our functions use spherical coordinates. The solution is use CoordinateTransform. As an example,

CoordinateTransform["Cartesian" -> "Spherical", {x, y, z}]

(*   {Sqrt[x^2 + y^2 + z^2], ArcTan[z, Sqrt[x^2 + y^2]], ArcTan[x, y]}  *)

We recognize this result as $(r,\theta,\phi)$. To make things as simple as possible, let's define our functions like this:

Clear[ele, f, g]

f[r_, θ_, ϕ_][r1_, ϕ1_] := 
    {r*Sin[θ]*Cos[ϕ] - r1*Cos[ϕ1], 
     r*Sin[θ]*Sin[ϕ] - r1*Sin[ϕ1],
     r*Cos[θ]}

g[r_, θ_, ϕ_][r1_, ϕ1_] := 
    Sqrt[(r*Sin[θ]*Cos[ϕ] - r1*Cos[ϕ1])^2 + 
          (r*Sin[θ]*Sin[ϕ] - r1*Sin[ϕ1])^2 +
          (r*Cos[θ])^2
    ]

ele[r_?NumericQ, θ_?NumericQ, ϕ_?NumericQ][R_?NumericQ] :=
    NIntegrate[(f[r, θ, ϕ][r1, ϕ1]/g[r, θ, ϕ][r1, ϕ1]^3)*r1,
      {ϕ1, 0, 2*π}, {r1, 0, R}, AccuracyGoal -> 8]

The reason we use those strange definitions is for compatibility with CoordinateTransformation. Now we can pass $(x,y,z)$ arguments to CoordinateTransformation, get back $(r,\theta,\phi)$, pass those to ele along with $R$ as a separate argument list and get back a Cartesian vector. Here is one way to plot the vectors:

With[{σ = 200, k = 9*10^9, R = 0.5},
 VectorPlot3D[
  (ele @@ CoordinateTransform["Cartesian" -> "Spherical", {x, y, z}])[R],
  {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]
 ]

What's going on? First, the @@ tells MMA to replace the head of List[r,θ, ϕ] with ele, so we end up with ele[r, θ, ϕ]. Then we follow that with [R] separately, all of which evaluates to our vector. That's the advantage of those strange definitions -- it lets us use @@ on the first 3 arguments only.

Also, we did not include $\sigma$ or $k$ in our functions and didn't use them in our With. That's because VectorPlot3D will rescale the vectors anyway, so it doesn't make any difference. So, really, the With is not necessary. It's only providing a value to $R$.