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][\[Phi]_, r_, \[Theta]_, \[Phi]1_, r1_] := {r*Sin[\[Theta]]*Cos[\[Phi]] - r1*Cos[\[Phi]1], r*Sin[\[Theta]]*Sin[\[Phi]] - r1*Sin[\[Phi]1], r*Cos[\[Theta]]}; \[ScriptR]Norm[\[Phi]_, r_, \[Theta]_, \[Phi]1_, r1_] := Sqrt[(r*Sin[\[Theta]]*Cos[\[Phi]] - r1*Cos[\[Phi]1])^2 + (r*Sin[\[Theta]]*Sin[\[Phi]] - r1*Sin[\[Phi]1])^2 + (r*Cos[\[Theta]])^2] // Simplify; ele[\[Phi]_?NumericQ, r_?NumericQ, \[Theta]_?NumericQ] := \[Sigma]*k* NIntegrate[(\[ScriptR][\[Phi], r, \[Theta], \[Phi]1, r1]/\[ScriptR]Norm[\[Phi], r, \[Theta], \[Phi]1, r1]^3)* r1, {\[Phi]1, 0, 2*\[Pi]}, {r1, 0, R}]; VectorPlot3D[ ele[\[Phi], r, \[Theta]] /. {\[Sigma] -> 200, k -> 9*10^9, R -> 0.5, r*Sin[\[Theta]]*Cos[\[Phi]] -> x, r*Sin[\[Theta]]*Sin[\[Phi]] -> y, r*Cos[\[Theta]] -> 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