# Dipolar magnetic field lines inside a cylinder

I'm drawing the magnetic field lines of a rotating dipole (magnetic field distorded by emission of radiation). Currently, my Mathematica 7.0 code is fully working, but I'm having a constraint which I don't know how to solve in a proper way: All the field lines I want to show should be tangent to the "light cylinder" surrounding the rotating dipole.

That means that the magnetic field vector $$\bf B$$ should be orthogonal to the cylinder's normal, located at $$x^2 + y^2 = 1$$ (the "light cylinder" have a radius of 1, in the code below) :

$$\quad \quad {\bf B \cdot n} = 0$$.

In the code below (for a simple non-rotating tilted dipole, to save space here), the field lines are starting on the surface of a sphere of radius omega < 1. The dipolar magnetic axis is tilted with an angle alpha, relative to the cylinder's axis. Take note that the constants alpha and omega are the two parameters of the field.

Using spherical coordinates theta and phi on the sphere (around the magnetic axis), I defined 30 starting points to grow the field lines from the sphere. The angle phi is uniformly distributed around the magnetic axis, but the theta angle is actually unknown. In the code below, I defined the theta angle by trial and error and this is what I want to fix.

So here's the code :

alpha=40Pi/180;

omega=2Pi (10000)/(299792458*0.001);

Mu0={Sin[alpha],0,Cos[alpha]};

r[x_,y_,z_]:=Sqrt[x^2+y^2+z^2]

Bdipolaire[x_,y_,z_]:=3(Mu0.{x,y,z}){x,y,z}/r[x,y,z]^5-Mu0/r[x,y,z]^3

NormeB[x_,y_,z_]:=Sqrt[Bdipolaire[x,y,z].Bdipolaire[x,y,z]]

Bx[x_,y_,z_]:={1,0,0}.Bdipolaire[x,y,z]
By[x_,y_,z_]:={0,1,0}.Bdipolaire[x,y,z]
Bz[x_,y_,z_]:={0,0,1}.Bdipolaire[x,y,z]

Angles:=
{0.1972815, 0.1988315,0.2014215,0.2049765,0.2090015,0.2127115,0.2151015,
0.2151015,0.2116515,0.2031515,0.1692915,0.1687315,0.1705015,0.1723815,
0.1742015,0.1759465,0.1776115,0.1792415,0.1808415,0.1824015,0.1839215,
0.1854085,0.1868585,0.1882765,0.1897015,0.1911415,0.1926035,0.1940415,
0.1953345,0.1963365}

NCourbes:=Length[Angles]

theta[n_]:= Angles[[n]]

phi[n_]:=(n-1)2Pi/NCourbes

CourbeMagnetique[n_]:=NDSolve[{
x'[s]==Bx[x[s],y[s],z[s]]/NormeB[x[s],y[s],z[s]],
y'[s]==By[x[s],y[s],z[s]]/NormeB[x[s],y[s],z[s]],
z'[s]==Bz[x[s],y[s],z[s]]/NormeB[x[s],y[s],z[s]],

x==omega (Sin[theta[n]]Cos[phi[n]]Cos[alpha]+Cos[theta[n]]Sin[alpha]),
y==omega Sin[theta[n]]Sin[phi[n]],
z==omega (Cos[theta[n]]Cos[alpha]-Sin[theta[n]]Cos[phi[n]]Sin[alpha])
},
{x,y,z}, {s,-1/omega,1/omega},
Method->Automatic,MaxSteps->10000000,
StoppingTest->(Sqrt[x[s]^2+y[s]^2+z[s]^2]<omega)]

Do[CourbeMagnetique[n],{n,1,NCourbes}]

Smin[n_]:=(x/.CourbeMagnetique[n])[][][][]
Smax[n_]:=(x/.CourbeMagnetique[n])[][][][]

GraphicSize:=1.5

GrapheMagnetique[n_]:=
ParametricPlot3D[
Evaluate[{x[s],y[s],z[s]}/.CourbeMagnetique[n]],
{s,Smin[n],Smax[n]},
PlotStyle->{Directive[AbsoluteThickness],Blue},
MaxRecursion->7,PerformanceGoal->"Quality"]

CercleLumiere:=
ParametricPlot3D[{Cos[p],Sin[p],0},{p,0,2Pi},
PlotStyle->{Directive[Thick,RGBColor[0.40,0.70,0.40]]},
PerformanceGoal->"Quality"]

AxesCartesiens=
{Line[GraphicSize {{-1,0,0},{1,0,0}}],
Line[GraphicSize {{0,-1,0},{0,1,0}}],
Line[GraphicSize {{0,0,-1},{0,0,1}}]};

AxeMagnetique=
Line[(2/3){{-Sin[alpha],0,-Cos[alpha]},{Sin[alpha],0,Cos[alpha]}}];

AxeRotation=Line[(1/3){{0,0,-1},{0,0,1}}];

AxesReference:=
Graphics3D[{
{Thin,GrayLevel[0.7],Dashed,AxesCartesiens},
{Thick,Blue,AxeMagnetique},
{Thick,RGBColor[0.40,0.70,0.40],AxeRotation}}]

Pulsar:=Graphics3D[Sphere[{0,0,0}, omega]]

Graphique=
Show[
Table[GrapheMagnetique[n],{n,1,NCourbes}],CercleLumiere,AxesReference,Pulsar,
PlotRange->{
{-GraphicSize,GraphicSize},
{-GraphicSize,GraphicSize},
{-GraphicSize,GraphicSize}},
Boxed->False,Axes->False,Lighting->"Neutral",
SphericalRegion->True,ViewPoint->{0,0,1}]


I would like to post a nice picture here to show what I'm trying to achieve, but that new picture upload interface don't allow me to upload a picture (what the hell!? Safari isn't supported?)

EDIT : Here's the picture I made with my handmade (madness) theta angles : And here's another one for another alpha angle: Pictures decription: The green circle is the light cylinder's equator. The green vertical axis is the symetry (rotation) axis. The blue line is the magnetic axis. Each curve on these pictures is tangent to the light cylinder. To achieve this, I need some proper theta angles for the field line starting point, around the magnetic axis.

Now, my question is this: how can I modify the code above so that Mathematica 7.0 calculates the proper 30 theta starting angles that makes the field lines tangent to the light cylinder, from the alpha and omega inputs?

Currently, I need to calculate all the theta angles myself using a painfull trial and error method, and was able to do it (very laboriously) for 5 alpha values (30, 40, 45, 60, 90 degrees) and 1 omega value only. This is madness!

EDIT 2 : For what it's worth, it can be shown that the theta angle that gives a field line tangent to the light cylinder is relatively close to $$\arcsin{\sqrt{\omega}}$$ :

$$\quad\quad \theta = \arcsin{\sqrt{\omega}} + \delta\theta$$.

When I've found the $$\theta$$ angles with my brute and painfull method, I was varying the angle around $$\arcsin{\sqrt{\omega}}$$ (for each value of $$\phi$$) and was evaluating the value of the scalar product $$\mathcal{P} \equiv \bf B \cdot \bf n$$ :

$$\quad\quad \mathcal{P} = B_x(\cos{\phi'}, \, \sin{\phi'}, \, z) \cos{\phi'} + B_y(\cos{\phi'}, \, \sin{\phi'}, z) \sin{\phi'}$$.

Here, $$\phi'$$ is the azimutal angle defined on the equator (or on the light cylinder), not around the magnetic axis, while $$\phi$$ is defined around the magnetic axis.

When $$\mathcal{P}$$ was far away from 0, I made another variation of $$\theta$$, until $$\mathcal{P}$$ was "close enough" to 0. Then I got the $$\theta$$ angle, for a given $$\phi$$. Remember that the two angles defined in my code above (theta $$\equiv \theta$$ and phi $$\equiv \phi$$) are defined on the sphere of radius omega $$\equiv \omega$$, around the tilted magnetic axis.

Maybe something similar could be done automatically with Mathematica ?

• "that new picture upload interface don't allow me to upload a picture" - you can try uploading to imgur.com's site, link to it here, and then somebody can edit in the picture for you. – J. M.'s torpor Aug 3 '15 at 14:47
• Ok, I'm trying this trick. But the older upload interface was working great before. Why the change ? – Cham Aug 3 '15 at 14:50
• I don't know either. FWIW, I detest the drag-and-drop business myself. – J. M.'s torpor Aug 3 '15 at 14:51
• I just pasted another picture, for another tilted angle. – Cham Aug 3 '15 at 14:56
• I am not having any problems with the new drag-and-drop image loader. I am running Safari 8.0.3 on Yosemite (10.10.2) – m_goldberg Aug 3 '15 at 15:25

## 1 Answer

The main issue is simply that your constraint should not be imposed after the integration of the field lines, but beforehand. This means that we should choose the starting points from which the differential equations of the field lines are integrated to lie on the desired cylinder right from the beginning. Then, all you have to do is to impose the tangentiality condition to points on the cylinder surface.

To illustrate things, I decided to re-use my own field line plotting functions from this answer, but I added the option to terminate field lines based on a WhenEvent condition. Therefore, I'll first post the new version of the general plotting code here (I'll add it to the earlier answer later):

fieldSolve::usage =
"fieldSolve[f,x,x0,\!$$\*SubscriptBox[\(t$$, $$max$$]\)] \
symbolically takes a vector field f with respect to the vector \
variable x, and then finds a vector curve r[t] starting at the point \
x0 satisfying the equation dr/dt=α f[r[t]] for \
t=0...\!$$\*SubscriptBox[\(t$$, $$max$$]\). Here α=1/|f[r[t]]| \
for normalization.";

fieldLinePlot[field_, varList_, seedList_, opts : OptionsPattern[]] :=
Module[{sols, localVars, var, localField, plotOptions, tubeFunction,
tubePlotStyle, postProcess = {}, stopEvent},
plotOptions = FilterRules[{opts}, Options[ParametricPlot3D]];
tubeFunction = OptionValue["TubeFunction"];
If[tubeFunction =!= None,
tubePlotStyle = Cases[OptionValue[PlotStyle], Except[_Tube]];
plotOptions =
FilterRules[plotOptions,
Except[{PlotStyle, ColorFunction, ColorFunctionScaling}]];
postProcess =
Line[x_] :>
Join[tubePlotStyle, {CapForm["Butt"],
Tube[x, tubeFunction @@@ x]}]];
If[Length[seedList[[1, 1]]] != Length[varList],
Print["Number of variables must equal number of initial conditions\
\nUSAGE:\n" <> fieldLinePlot::usage]; Abort[]];
localVars = Array[var, Length[varList]];
stopEvent =
"StoppingEvent" /.
ReleaseHold[
FilterRules[Unevaluated[{opts}], "StoppingEvent"] /.
Thread[Map[HoldPattern, Unevaluated[varList]] ->
localVars]] /. Options[fieldLinePlot, "StoppingEvent"];
localField =
ReleaseHold[
Hold[field] /.
Thread[Map[HoldPattern, Unevaluated[varList]] -> localVars]];
(*Assume that each element of seedList specifies a point AND the \
length of the field line:*)
Show[ParallelTable[
ParametricPlot3D[
Evaluate[
Through[#[t]]], {t, #[[1, 1, 1, 1]], #[[1, 1, 1, 2]]},
Evaluate@Apply[Sequence, plotOptions]] &[
fieldSolve[localField, localVars, seedList[[i, 1]],
seedList[[i, 2]], "StoppingEvent" -> stopEvent]] /.
postProcess, {i, Length[seedList]}]]]

Options[fieldSolve] = {"StoppingEvent" -> False};

fieldSolve[field_, varlist_, xi0_, tmax_, OptionsPattern[]] :=
Module[{xiVec, equationSet, t, stopEvent},
If[Length[varlist] != Length[xi0],
Print["Number of variables must equal number of initial conditions\
\nUSAGE:\n" <> fieldSolve::usage]; Abort[]];
xiVec = Through[varlist[t]];
(*Below,Simplify[
equationSet] would cost extra time and doesn't help with the \
numerical solution,so don't try to simplify.*)
stopEvent = OptionValue["StoppingEvent"] /. Thread[varlist -> xiVec];
equationSet =
Join[Thread[
Map[D[#, t] &, xiVec] ==
Normalize[field /. Thread[varlist -> xiVec]]],
Thread[(xiVec /. t -> 0) == xi0]];
(*This is where the differential equation is solved.The Quiet[] \
command suppresses warning messages because numerical precision isn't \
crucial for our plotting purposes:*)
Map[Head,
First[xiVec /.
Quiet[NDSolve[
Append[equationSet,
WhenEvent[Evaluate@stopEvent, "StopIntegration"]],
xiVec, {t, 0, tmax}]]], 2]]

Options[fieldLinePlot] =
Join[Options[ParametricPlot3D], {"TubeFunction" -> None,
"StoppingEvent" -> False}];

SyntaxInformation[
fieldLinePlot] = {"LocalVariables" -> {"Solve", {2, 2}},
"ArgumentsPattern" -> {_, _, _, OptionsPattern[]}};

SetAttributes[fieldSolve, HoldAll];
SetAttributes[fieldLinePlot, HoldAll];


With this, we can make the plots for your application as follows. I changed the order of the definitions so that alpha remains as a symbolic parameter in all the functions, and the list of starting points. The latter is obtained by using Solve[tangentCondition==0,z] on a point on the cylinder, whose coordinates x and y are determined by an angle ϕ. There are two solutions for z, and I use both. Also, I go in both possible directions from each point obtained in this way. In my function fieldLinePlot, each element of the point list is expected to come with a maximum length up to which NDSolve will trace the field line. These lengths are set to 10 for all starting points. The direction in which to integrate is specified by the sign of these lengths.

Clear[alpha, omega, Mu0]

omega = 2 Pi (10000)/(299792458*0.001);

Mu0 = {Sin[alpha], 0, Cos[alpha]};

r[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2]

Bdipolaire[x_, y_, z_] :=
3 (Mu0.{x, y, z}) {x, y, z}/r[x, y, z]^5 - Mu0/r[x, y, z]^3

tangentCondition =
FullSimplify[(Bdipolaire[x, y, z].{Cos[ϕ], Sin[ϕ],
0}) /. {x -> Cos[ϕ], y -> Sin[ϕ]}, ϕ > 0];

{zStart1, zStart2} = z /. Solve[tangentCondition == 0, z];

AxeMagnetique =
Tube[Line[(2/3) {{-Sin[alpha], 0, -Cos[alpha]}, {Sin[alpha], 0,
Cos[alpha]}}], .05];

pointList[alphaVar_] :=
Module[{length = 10},
N[Flatten[
Table[{{{Cos[ϕ], Sin[ϕ], zStart1},
length}, {{Cos[ϕ], Sin[ϕ], zStart2},
length}, {{Cos[ϕ], Sin[ϕ],
zStart1}, -length}, {{Cos[ϕ], Sin[ϕ],
zStart2}, -length}}, {ϕ, Pi/5, 2 Pi, Pi/5}],
1]] /. alpha -> alphaVar]

Block[{alpha = 1.1},
Show[
fieldLinePlot[
Bdipolaire[x, y, z], {x, y, z}, pointList[alpha],
PlotStyle -> {Cyan, Specularity[White, 16]}, PlotRange -> All,
Boxed -> False, Axes -> None,
"StoppingEvent" -> (x^2 + y^2 + z^2 < omega^2 || x^2 + y^2 > 1)],
Graphics3D[{Orange, AxeMagnetique, White,
Sphere[{0, 0, 0}, omega]}],
Graphics3D[{Opacity[.5], White,
Cylinder[{{0, 0, -2}, {0, 0, 2}}, 1]}], Background -> Black,
PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}, Lighting -> "Neutral"]] Block[{alpha = Pi/2},
Show[
fieldLinePlot[
Bdipolaire[x, y, z], {x, y, z}, pointList[alpha],
PlotStyle -> {Cyan, Specularity[White, 16]}, PlotRange -> All,
Boxed -> False, Axes -> None,
"StoppingEvent" -> (x^2 + y^2 + z^2 < omega^2 || x^2 + y^2 > 1)],
Graphics3D[{Orange, AxeMagnetique, White,
Sphere[{0, 0, 0}, omega]}],
Graphics3D[{Opacity[.5], White,
Cylinder[{{0, 0, -2}, {0, 0, 2}}, 1]}], Background -> Black,
PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}, Lighting -> "Neutral"]] Fix for version 7:

Because WhenEvent doesn't exist in versions 7 and 8, you could instead use the following version of fieldSolve. It uses the StoppingTest of NDSolve instead, as in the question.

fieldSolve[field_, varlist_, xi0_, tmax_, OptionsPattern[]] :=
Module[{xiVec, equationSet, t, stopEvent},
If[Length[varlist] != Length[xi0],
Print["Number of variables must equal number of initial conditions\
\nUSAGE:\n" <> fieldSolve::usage]; Abort[]];
xiVec = Through[varlist[t]];
(*Below,Simplify[
equationSet] would cost extra time and doesn't help with the \
numerical solution,so don't try to simplify.*)
stopEvent = OptionValue["StoppingEvent"] /. Thread[varlist -> xiVec];
equationSet =
Join[Thread[
Map[D[#, t] &, xiVec] ==
Normalize[field /. Thread[varlist -> xiVec]]],
Thread[(xiVec /. t -> 0) == xi0]];
(*This is where the differential equation is solved.The Quiet[] \
command suppresses warning messages because numerical precision isn't \
crucial for our plotting purposes:*)
Map[Head,
First[xiVec /.
Quiet[NDSolve[equationSet, xiVec, {t, 0, tmax},
StoppingTest -> stopEvent]]], 2]]


Furthermore, in Mathematica version 7 the stopping conditions are apparently not recognized properly when the integration goes beyond the condition for being inside the cylinder. As a workaround, you just have to be more judicious in the choice of initial conditions for the field lines, making sure that you pick only the ones that are guaranteed to be inside the cylinder, not just being tangential to it. There are tangent field lines with a large radius of curvature that "hug" the cylinder from the outside, and they have to be excluded. Because of their large radius, these correspond to large z values. I find that they can be eliminated by changing the list of starting points to

pointList[alphaVar_] :=
Module[{length = 10},
N[Flatten[
Table[{{{Cos[ϕ], Sin[ϕ], zStart1},
length}, {{Cos[ϕ], Sin[ϕ],
zStart1}, -length}}, {ϕ, Pi/5, 2 Pi, Pi/5}], 1]] /.
alpha -> alphaVar]


in the code listed above. This uses only the first solution of the tangency condition, zStart1. Dropping zStart2 also eliminates some legitimate field lines that are entirely within the cylinder, but it leaves enough field lines to give a nice plot.

Modifying the original code of the question to use the above points

The original code isn't modularized and causes many unnecessary evaluations because it uses := inappropriately and does a futile Do loop. That's why I didn't use it initially. I fixed those problems, and added my initial conditions. Here is the result that I think is closest to the code in the question:

Clear[alpha, omega, Mu0]

omega = 2 Pi (10000)/(299792458*0.001);

Mu0 = {Sin[alpha], 0, Cos[alpha]};

r[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2]

Bdipolaire[x_, y_, z_] :=
3 (Mu0.{x, y, z}) {x, y, z}/r[x, y, z]^5 - Mu0/r[x, y, z]^3

tangentCondition =
FullSimplify[(Bdipolaire[x, y, z].{Cos[ϕ], Sin[ϕ],
0}) /. {x -> Cos[ϕ], y -> Sin[ϕ]}, ϕ > 0];

{zStart1, zStart2} = z /. Solve[tangentCondition == 0, z];

alpha = 40 Pi/180;

omega = 2 Pi (10000)/(299792458*0.001);

NormeB[x_, y_, z_] := Sqrt[Bdipolaire[x, y, z].Bdipolaire[x, y, z]]

Bx[x_, y_, z_] := {1, 0, 0}.Bdipolaire[x, y, z]
By[x_, y_, z_] := {0, 1, 0}.Bdipolaire[x, y, z]
Bz[x_, y_, z_] := {0, 0, 1}.Bdipolaire[x, y, z]

NCourbes = 5;

phi[n_] := (n - 1) 2 Pi/NCourbes

Table[
CourbeMagnetique[n] =
NDSolve[{x'[s] == Bx[x[s], y[s], z[s]]/NormeB[x[s], y[s], z[s]],
y'[s] == By[x[s], y[s], z[s]]/NormeB[x[s], y[s], z[s]],
z'[s] == Bz[x[s], y[s], z[s]]/NormeB[x[s], y[s], z[s]],
x == Cos[phi[n]],
y == Sin[phi[n]],
z == zStart1 /. ϕ -> phi[n]},
{x, y, z}, {s, -10, 10}, Method -> Automatic,
MaxSteps -> 10000000,
StoppingTest -> (Sqrt[x[s]^2 + y[s]^2 + z[s]^2] < omega)],
{n, 1, NCourbes}];

Smin[n_] := (x /. CourbeMagnetique[n])[][][][]
Smax[n_] := (x /. CourbeMagnetique[n])[][][][]

GraphicSize := 1.5

GrapheMagnetique[n_] :=
ParametricPlot3D[
Evaluate[{x[s], y[s], z[s]} /. CourbeMagnetique[n]], {s, Smin[n],
Smax[n]}, PlotStyle -> {Directive[AbsoluteThickness], Blue},
MaxRecursion -> 7, PerformanceGoal -> "Quality"]

CercleLumiere =
ParametricPlot3D[{Cos[p], Sin[p], 0}, {p, 0, 2 Pi},
PlotStyle -> {Directive[Thick, RGBColor[0.40, 0.70, 0.40]]},
PerformanceGoal -> "Quality"];

AxesCartesiens = {Line[GraphicSize {{-1, 0, 0}, {1, 0, 0}}],
Line[GraphicSize {{0, -1, 0}, {0, 1, 0}}],
Line[GraphicSize {{0, 0, -1}, {0, 0, 1}}]};

AxeMagnetique =
Line[(2/3) {{-Sin[alpha], 0, -Cos[alpha]}, {Sin[alpha], 0,
Cos[alpha]}}];

AxeRotation = Line[(1/3) {{0, 0, -1}, {0, 0, 1}}];

AxesReference =
Graphics3D[{{Thin, GrayLevel[0.7], Dashed, AxesCartesiens}, {Thick,
Blue, AxeMagnetique}, {Thick, RGBColor[0.40, 0.70, 0.40],
AxeRotation}}];

Pulsar = Graphics3D[Sphere[{0, 0, 0}, omega]];

Graphique =
Show[Table[GrapheMagnetique[n], {n, 1, NCourbes}], CercleLumiere,
AxesReference, Pulsar,
PlotRange -> {{-GraphicSize, GraphicSize}, {-GraphicSize,
GraphicSize}, {-GraphicSize, GraphicSize}}, Boxed -> False,
Axes -> False, Lighting -> "Neutral", SphericalRegion -> True,
ViewPoint -> {0, 0, 1}]

• Jens, I tried your code above, and it doesn't compile well. I'm getting tons of error messages and the graphics doesn't combine. The first bloc of code compiles, apparently (there's no output, though). The second bloc gives the error messages. Take note that I'm working with Mathematica 7.0, so maybe there's a command that is not recognised by this version ? – Cham Aug 4 '15 at 0:33
• @Cham Right, the WhenEvent command doesn't exist in version 7. I can change it - but won't have time until a few hours from now. Maybe you can adapt your own code - the central ingredients are tangentCondition and zStart1, zStart2. – Jens Aug 4 '15 at 0:36
• Well, to be honest, I don't understand much of your code. I'm not a novice user of Mathematica, but I'm always at the "basics" level. I can wait for your final answer. There's no rush. – Cham Aug 4 '15 at 0:41
• I updated the answer with a fix for the older versions of Mathematica - hopefully this will work (I can't test on version 7). – Jens Aug 4 '15 at 5:17
• Is "EventLocator" no longer available in version 10? – J. M.'s torpor Aug 4 '15 at 6:27