Spatial rotation induces an unphysical blow up of the B field at the centre of a Halbach octupole

[Sylk1-9]

Magpylib version

4.x (Latest)

What happened?

Issue :

When computing the magnetic B field right in the centre of a Halbach octupole that has been spatially rotated, the result blows up.
The expected magnetics field should be zero though.

Code example :

The octupole is made of 32 individual CylinderSegment, distributed circularly.
The following code example build such octupole, with z as the axis of symmetry.
The B field vector is then computed for values distributed along the z axis. The result is close to zero, as expected.

However, rotating the same octupole around the y axis, then computing the B field vector for values distributed along the x axis (the new axis of symmetry), and the B field amplitude blows up to unphysical values.

Precisions and tests :

• Using the magpy routine rotate_from_angax or the scipyspatial.transform.Rotation routine from_rotvec provide similar results

• Strangely, playing with the radii, height, or magnetization of the disc segments doesn't change the position at which the B field explodes. It is always located around 50mm before or after the centre of the octupole, along the axis of symetry.

• The issue only appears when computing the magnetic field on the symmetry axis. Outside, the results between both orientations are similar.

Code example

import magpylib as magpy
import numpy as np
from scipy.spatial.transform import Rotation as R
from matplotlib import pyplot as plt

def buildOctupole():
    N = 32 # number of Segments
    phi = 360/N # angle between semgents
    r1 = 1
    r2 = 2
    z = 1 # heigh of segment
    M = 1 # Magnetization
    octupole = magpy.Collection()
    for deg in np.arange(0, 360, phi):
        rad = np.deg2rad(deg)
        magnetization = (M*np.cos(5*rad), M*np.sin(5*rad), 0)
        dimension = (r1, r2, z, deg-phi/2, deg+phi/2)
        position = (0, 0, 0)
        orientation = None # R.from_rotvec((0, 90, 0), degrees=True)
        color = 'rgb(%i, %i, %i)' % tuple(np.random.randint(0, 255, 3))
        segment = magpy.magnet.CylinderSegment(magnetization, dimension, position, orientation=orientation, style={'color':color})
        octupole.add(segment)
    return octupole



octupole = buildOctupole()

zline = np.zeros((100, 3))
zline[:, 2] = np.linspace(-100, 100, zline.shape[0])
B_zline = octupole.getB(zline)

xline = np.zeros((100, 3))
xline[:, 0] = np.linspace(-100, 100, xline.shape[0])
B_xline = octupole.rotate_from_angax(90, axis='y',anchor=0).getB(xline)


octupole.show()
plt.plot(zline[:, 2], np.abs(B_zline), label=("|Bx| z axis", "|By| z axis", "|Bz| z axis"))
plt.plot(xline[:, 0], np.abs(B_xline), label=("|Bx| x axis", "|By| x axis", "|Bz| x axis"), ls="--")

plt.yscale("log")
plt.ylabel("log(|B|)")
plt.xlabel("x | z")
plt.legend()
plt.show()

Additional context

• Ubuntu 22,
• python 3.10.6
• magpy 4.2.0
• numpy 1.24.2
• scipy 1.8.1

[Alexboiboi]

Hi @Sylk1-9,

Thanks for your interest in the library and for taking the time to explain the issue welll ;)

Looks like there is a numerical instability in the CylinderSegment computation.

I can reproduce it with a more minimal example:

import magpylib as magpy
import numpy as np
from matplotlib import pyplot as plt

cyl_seg = magpy.magnet.CylinderSegment(
    magnetization=[1000.0, 0.0, 0.0], dimension=[1.0, 2.0, 1.0, -10, 10]
)
zline = np.linspace((0, 0, -100), (0, 0, 100), 100)
B_zline = cyl_seg.getB(zline)

cyl_seg.rotate_from_angax(90, axis="y", anchor=0)
xline = zline[:, ::-1]
B_xline = cyl_seg.getB(xline)

plt.plot(
    zline[:, 2],
    np.abs(B_zline),
    label=("|Bx| z axis", "|By| z axis", "|Bz| z axis"),
)
plt.plot(
    zline[:, 2],
    np.abs(B_xline),
    label=("|Bx| x axis", "|By| x axis", "|Bz| x axis"),
    ls="--",
)

plt.yscale("log")
plt.ylabel("log(|B|)")
plt.xlabel("x | z")
plt.legend()
plt.show()

If the line is a sensor which is rotated, there is no issue though:

import magpylib as magpy
import numpy as np
from matplotlib import pyplot as plt

cyl_seg = magpy.magnet.CylinderSegment(
    magnetization=[1000.0, 0.0, 0.0], dimension=[1.0, 2.0, 1.0, -10, 10]
)
pixel = np.linspace((0, 0, -100), (0, 0, 100), 100)
sens = magpy.Sensor(pixel=pixel)
B_zline = cyl_seg.getB(sens)

cyl_seg.rotate_from_angax(90, axis="y", anchor=0)
sens.rotate_from_angax(90, axis="y", anchor=0)
B_xline = cyl_seg.getB(sens)

plt.plot(
    pixel[:, 2],
    np.abs(B_zline),
    label=("|Bx| z axis", "|By| z axis", "|Bz| z axis"),
)
plt.plot(
    pixel[:, 2],
    np.abs(B_xline),
    label=("|Bx| x axis", "|By| x axis", "|Bz| x axis"),
    ls="--",
)

plt.yscale("log")
plt.ylabel("log(|B|)")
plt.xlabel("x | z")
plt.legend()
plt.show()

@OrtnerMichael, can you have take at this?

[OrtnerMichael]

Hi @Sylk1-9

Not sure what is meant by "blowing up". When I run your code I get this:

What do you get ?

@Alexboiboi

When I run your first code I get this

and not what you show. Your second code does not compile - its missing the zline. But when I add

zline = np.linspace((0, 0, -100), (0, 0, 100), 100)

i get the same as you do. You sure about the first code? What operating system did you run this on ?

[Sylk1-9]

@OrtnerMichael

Hum, indeed, I should have included the plot directly. Mine is very different from yours, when I run the code of my initial post.

About @Alexboiboi code examples, I get the exact two same figures as he got as well.

[Alexboiboi]

zline = np.linspace((0, 0, -100), (0, 0, 100), 100)

forgot to remove the line. Fixed it in the comment now.
I tried my codes on Mac OS and Windows with the same result

[Alexboiboi]

@OrtnerMichael

Seems to be an issue of recognizing the right CylinderSegment cases.

magpylib/magpylib/_src/fields/field_BH_cylinder_segment.py

Lines 2416 to 2421 in 6cd0211

	# transform obs_pos to Cy CS --------------------------------------------
	x, y, z = observers.T
	r, phi = np.sqrt(x**2 + y**2), np.arctan2(y, x)
	pos_obs_cy = np.concatenate(((r,), (phi,), (z,)), axis=0).T
	# determine when points lie inside and on surface of magnet -------------

adding

observers[np.abs(observers)<1e-12] = 0.

between L2417-L2418 solves the issue. It may need deeper analysis for robustness though...

[OrtnerMichael]

agreed !!! lets keep this issue open until this is fixed !