The following stripped down code works, but I'm experiencing some strange warnings: from an old document made with Mathematica 7.0 (but used in version 13.2), I'm getting several warning messages about "Event location failed to converge to requested accuracy or precision...". The output is okay, though. When I copy paste the exact same code in a new document made with Mathematica 13.2, I don't get any warning/error message, but the StoppingTest is drawn in red. So the issue lies there (indicated with the comment (* ISSUE IS HERE! *) in the code below).
MomMagn = {1, 1, 1}; (* norm of each magnetic moment *) thetaMu = {45, 45, 0}Pi/180; (* inclinaison relative to Z axis *) phiMu = {90, 270, 0}Pi/180; (* azimutal angle around Z axis *) source[k_] := {{3, 0, 0}, {-8, 0, 0}, {0, 8, 0}}[[k]] (* position of each magnetic sphere *) rotationX[k_] := {{1, 0, 0}, {0, Cos[thetaMu[[k]]], -Sin[thetaMu[[k]]]}, {0, Sin[thetaMu[[k]]], Cos[thetaMu[[k]]]}} (* rotation matrix around X axis *) rotationZ[k_] := {{Cos[phiMu[[k]]], -Sin[phiMu[[k]]], 0}, {Sin[phiMu[[k]]], Cos[phiMu[[k]]], 0}, {0, 0, 1}} (* rotation matrix around Z axis *) mu[k_] := MomMagn[[k]]rotationZ[k].rotationX[k].{0, 0, 1} (* magnetic moments *) Nspheres = Length[MomMagn]; (* number of magnetic spheres *) (* magnetic field definition: *) xyz[s_] := {x[s], y[s], z[s]} r[s_, k_] := Norm[xyz[s] - source[k]] dipole[s_, k_] := 3 mu[k].(xyz[s] - source[k])(xyz[s] - source[k])/r[s, k]^5 - mu[k]/r[s, k]^3 magneticfield[s_] := Sum[dipole[s, k], {k, 1, Nspheres}] (* magnetic field lines definition: *) Nlines = 15; theta[n_] := Which[n == -2, 20, n == -1, 0, n == 0, 0, n == 1, 20, n == 2, 27.5]Pi/180 phi[n_, m_] := If[n > 0, m 2 Pi/Nlines, 0] sphericalstart[n_, m_] := If[n < 0, -1, 1] {Sin[theta[n]] Cos[phi[n, m]], Sin[theta[n]] Sin[phi[n, m]], Cos[theta[n]]} linestart[n_, m_, k_] := rotationZ[k].rotationX[k].sphericalstart[n, m] + source[k] magneticCurves[n_, m_, k_] := NDSolve[{ x'[s] == If[n < 0, -1, 1]{1, 0, 0}.magneticfield[s]/Norm[magneticfield[s]], y'[s] == If[n < 0, -1, 1]{0, 1, 0}.magneticfield[s]/Norm[magneticfield[s]], z'[s] == If[n < 0, -1, 1]{0, 0, 1}.magneticfield[s]/Norm[magneticfield[s]], x[0] == {1, 0, 0}.linestart[n, m, k], y[0] == {0, 1, 0}.linestart[n, m, k], z[0] == {0, 0, 1}.linestart[n, m, k] }, {x, y, z}, {s, 0, 500}, Method -> Automatic, MaxSteps -> Automatic, StoppingTest -> (Apply[Or, Table[EuclideanDistance[xyz[s], source[kk]] < 1, {kk, 1, Nspheres}]] || r[s, k] > 25)] (* ISSUE IS HERE! *) Do[magneticCurves[n, m, k], {k, 1, Nspheres}, {n, -2, 2}, {m, 1, Nlines}] Smin[n_, m_, k_] := (x/.magneticCurves[n, m, k])[[1]][[1]][[1]][[1]] (* start of each field line *) Smax[n_, m_, k_] := (x/.magneticCurves[n, m, k])[[1]][[1]][[1]][[2]] (* end of each field line *) magneticGraph[n_, m_, k_] := ParametricPlot3D[Evaluate[xyz[s]/.magneticCurves[n, m, k]], {s, Smin[n, m, k], Smax[n, m, k]}] Show[Table[magneticGraph[n, m, k], {n, -2, 2}, {m, 1, Nlines}, {k, 1, Nspheres}], PlotRange -> All] Preview of what the code is doing:
So what should I use instead of StoppingTest? If I remove completely this part, the output isn't right at all. The field lines output is right when I still use the red StoppingTest line. I don't understand what's going on here.
Note: This code is slow (on a Silicon Mac M2 Pro with 32GB ram). This is NOT a concern to me. The issue isn't about the compilation speed and I don't need to optimize the whole code.
