0
$\begingroup$

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

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. enter image description here

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!

1 https://arxiv.org/pdf/0909.3880.pdf

\[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}]
Improve this question
$\endgroup$
3
  • $\begingroup$ Do your z+ and z- need to be squared inside the square root in definitions of alpha and beta? $\endgroup$
    – Rudy Potter
    Commented Oct 25, 2022 at 15:52
  • $\begingroup$ Not sure why you got rid of the outer square root in k2+ and k2- $\endgroup$
    – Rudy Potter
    Commented Oct 25, 2022 at 16:16
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ariana Fenris
    Commented Oct 26, 2022 at 0:28

1 Answer 1

Reset to default
1
$\begingroup$

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.