# Empty VectorPlot3D in Magnetic Field of Solenoid

I have been trying to plot the 3D magnetic field of a solenoid, with the below expression 1:

With C, the generalised complete elliptic integral, defined as follows:

I have implemented all of this in the below code (I apologise for the length), in which I also converted to Cartesian coordinates to match the coordinate system of VectorPlot3D but when I run VectorPlot3D the plot is simply empty. In debugging I output TraditionalForm[NewBx[x,y,z]], which returns a strange conditional that I thought could be the cause of the empty plot.

The expressions are undefined at rho=r1 and z=+-length. I attempted to exclude these points by making newBx, newBy, newBz Piecewise functions but it does not seem to work either.

What am I doing wrong? Thank you very much!

\[Mu]0 = 4*Pi*10^-7*10^2;

(*coil properties*)

r1 = 3/2;(*radius*)
current = 2;
length = 1;

(*define helper functions*)

zplus[z_] := z + length;
zminus[z_] := z - length;
\[Alpha]plus[z_, \[Rho]_] := r1/Sqrt[zplus[z] + (\[Rho] + r1)^2];
\[Alpha]minus[z_, \[Rho]_] := r1/Sqrt[zminus[z] + (\[Rho] + r1)^2];
\[Beta]plus[z_, \[Rho]_] := zplus[z]/Sqrt[zplus[z] + (\[Rho] + r1)^2];
\[Beta]minus[z_, \[Rho]_] := zminus[z]/Sqrt[zminus[z] + (\[Rho] + r1)^2];
\[Gamma][\[Rho]_] := (r1 - \[Rho])/(r1 + \[Rho]);
k2plus[z_, \[Rho]_] := (zplus[z]^2 + (r1 - \[Rho])^2)/(zplus[z]^2 + (r1 + \[Rho])^2);(*NOTE: removing the Sqrt[] cause it always comes ^2*)

k2minus[z_, \[Rho]_] := (zminus[z]^2 + (r1 - \[Rho])^2)/(zminus[z]^2 + (r1 + \[Rho])^2);
Cellip[kc_, p_, c_, s_] := c*CarlsonRF[0, kc, 1] + 1/3*(s - p*c)*CarlsonRJ[0, kc, 1, p];(*complete general elliptic integral*)

B0 = (\[Mu]0*10*current)/Pi;

newBx[x_, y_, z_] := (B0*(\[Alpha]plus*Cellip[k2plus, 1, 1, -1] - \[Alpha]minus*Cellip[k2minus, 1, 1, -1]))*x/Sqrt[x^2 + y^2] /. \[Rho] -> Sqrt[x^2 + y^2]
newBy[x_, y_, z_] := (B0*(\[Alpha]plus*Cellip[k2plus, 1, 1, -1] - \[Alpha]minus*Cellip[k2minus, 1, 1, -1]))*y/Sqrt[x^2 + y^2] /. \[Rho] -> Sqrt[x^2 + y^2]
newBz[x_, y_, z_] := (B0*r1)/(r1 + \[Rho])*(\[Beta]plus*Cellip[k2plus, \[Gamma]^2, 1, \[Gamma]] - \[Beta]minus*Cellip[k2minus, \[Gamma]^2, 1, \[Gamma]]) /. \[Rho] -> Sqrt[x^2 + y^2]

VectorPlot3D[{Evaluate@newBx[x, y, z], Evaluate@newBy[x, y, z], Evaluate@newBz[x, y, z]}, {x, 0, 5}, {y, 0, 5}, {z, -length - 2, length + 2}]

• Do your z+ and z- need to be squared inside the square root in definitions of alpha and beta? Commented Oct 25, 2022 at 15:52
• Not sure why you got rid of the outer square root in k2+ and k2- Commented Oct 25, 2022 at 16:16
• Gosh @RudyPotter sorry about the error I fixed the z+ and z- they should be squared. But the plot remains empty. I removed the outer square root in k2+ and k2- because the input to the elliptic integrals is always k2^2. Commented Oct 26, 2022 at 0:28

## 1 Answer

I found the error, I had forgotten to enter the arguments of the various helper functions ie [z,\[Rho] etc. It plots beautifully now: Thank you so much everyone!