Magnetic field is spiking but with rather not random pattern.

[cwahn]

What is your question?

Hello guys, I'm quite new to the magpylib.

I am simulating the magnetic field of a radially magnetized ring-shaped magnet. I followed the example and I managed to find the magnetic field plot as I expected.

The design problem has strict conditions regarding the minimum and maximum magnetic field intensity. So I tried to visualize and also analyze the characteristics of magnetic fields on a certain path.

However, the problem is when I plot the magnetic field, it's spiking and rather not a random pattern when it seems smooth on visualization. Doesn't seem like a quantization error or an actual issue of the library, but I cannot find a reason.

What could be possibly wrong?

from typing import List, Tuple, Optional
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt
import numpy as np
from math import sqrt

import magpylib as magpy

# Global variables

br_mT = {
    "N35": 1210,
    "N45": 1350,
    "N52": 1450,
}

material = "N45"

br_mT = br_mT[material]

magnetization = np.array((br_mT, 0, 0))
inner_diameter_m = 32e-3
outer_diameter_m = 40e-3
height_m = 10e-3

min_B_value = 30  # minimum B field value in mT you want to display
max_B_value = 1500  # maximum B field value in mT you want to display

grid_n = 1000  # number of points in the grid

def boundary_squre_length(outer_diameter_m: float) -> float:
    return outer_diameter_m * 1.5

def point_index_from_position(
    boudnary_l: float, p: Tuple[float, float]
) -> Tuple[int, int]:
    x_idx = int((p[0] + boudnary_l / 2) / boudnary_l * grid_n)
    y_idx = int((p[1] + boudnary_l / 2) / boudnary_l * grid_n)
    return x_idx, y_idx

def circle_points(r: float, n: int) -> List[Tuple[float, float]]:
    points = map(
        lambda i: (
            r * np.cos(2 * np.pi * i / n),
            r * np.sin(2 * np.pi * i / n),
        ),
        range(n),
    )
    return list(points)

def circle_indices(r: float, n: int) -> List[Tuple[int, int]]:
    l = boundary_squre_length(r * 2)
    ps = circle_points(r, n)
    idxs = map(lambda p: point_index_from_position(l, p), ps)
    return list(idxs)

def circle_bs(b: np.ndarray, r: float, n: int) -> List[Tuple[float, float, float]]:
    idxs = circle_indices(r, n)
    bs = map(lambda idx: b[idx[0], idx[1]], idxs)
    b_values = map(lambda b: (b[0], b[1], b[2]), bs)
    return list(b_values)

def plot_path_b_xy(bs: List[Tuple[float, float, float]]):
    n = len(bs)
    degs = list(map(lambda i: i / n * 360, range(n)))

    bxs_mT = list(map(lambda b: b[0] * 1e3, bs))
    bys_mT = list(map(lambda b: b[1] * 1e3, bs))

    # ! temp
    print(f"bxs: {bxs_mT}")

    fig, ax = plt.subplots(figsize=(16, 9))
    ax.plot(degs, bxs_mT, label="Bx")
    ax.plot(degs, bys_mT, label="By")

    ax.legend()
    ax.set_title("B_x and B_y on path ")
    ax.set_xlabel("degrees")
    ax.set_ylabel("B (mT)")

    plt.show()

def max_b_rad_T(bs: List[Tuple[float, float, float]]) -> float:
    max_b = max(bs, key=lambda b: abs(b[1]))
    return max_b[1]

def max_b_tan_T(bs: List[Tuple[float, float, float]]) -> float:
    max_b = max(bs, key=lambda b: abs(b[0]))
    return max_b[0]

def min_b_xy_T(bs: List[Tuple[float, float, float]]) -> float:
    b_xy_abss = list(map(lambda b: sqrt(b[0] ** 2 + b[1] ** 2), bs))
    return min(b_xy_abss)

def plot_b_field(
    b: np.ndarray, x_grid: np.ndarray, y_grid: np.ndarray, r: Optional[float] = None
):
    # Create a Matplotlib figure
    fig, ax = plt.subplots(figsize=(16, 16))

    # levels = np.linspace(min_B_value, max_B_value, 100)
    levels = np.geomspace(min_B_value, max_B_value, 100)

    normB = np.linalg.norm(b, axis=2)

    # combine streamplot with contourf
    cp = ax.contourf(
        x_grid,
        y_grid,
        normB * 1000,  # T->mT
        cmap="rainbow",
        levels=levels,
        norm=LogNorm(vmin=min_B_value, vmax=max_B_value),
        zorder=1,
        extend="neither",  # Do not extend the colorbar to indicate out-of-range values
    )

    ax.streamplot(
        x_grid,
        y_grid,
        b[:, :, 0],
        b[:, :, 1],
        color="w",
        density=1.5,
        linewidth=1,
        zorder=3,
    )

    # Add colorbar
    fig.colorbar(cp, ax=ax, label="|B| (mT)", format="%0.1f")

    # Outline magnet boundary
    ts = np.linspace(0, 2 * np.pi, 50)
    ax.plot(
        outer_diameter_m / 2 * np.cos(ts),
        outer_diameter_m / 2 * np.sin(ts),
        "w-",
        lw=2,
        zorder=2,
    )
    ax.plot(
        inner_diameter_m / 2 * np.cos(ts),
        inner_diameter_m / 2 * np.sin(ts),
        "w-",
        lw=2,
        zorder=2,
    )

    if r:
        ax.plot(
            r * np.cos(ts),
            r * np.sin(ts),
            "k--",
            lw=2,
            zorder=2,
        )

    # Figure styling
    ax.set(
        xlabel="x-position (m)",
        ylabel="z-position (m)",
        aspect=1,
    )

    ax.set_title(
        f"Ring magnet with material {material}\n inner diameter {inner_diameter_m * 1000:.0f} mm\n outer diameter {outer_diameter_m * 1000:.0f} mm\n height {height_m * 1000:.0f} mm, in z=0 plane"
    )

    plt.tight_layout()
    plt.show()

 # Create an observer grid in the xy-symmetry plane
l = boundary_squre_length(outer_diameter_m)
xs = np.linspace(-l / 2, l / 2, grid_n)
ys = np.linspace(-l / 2, l / 2, grid_n)
x_grid, y_grid = np.meshgrid(xs, ys)
grid = np.stack([x_grid, y_grid, np.zeros((grid_n, grid_n))], axis=2)

# Compute magnetic field on grid - using the functional interface
b_field = magpy.getB(
    "CylinderSegment",
    observers=grid.reshape(-1, 3),
    dimension=(inner_diameter_m / 2, outer_diameter_m / 2, height_m / 2, 0, 360),
    polarization=(br_mT * 1e-3, 0, 0),
)

b_field = b_field.reshape(grid.shape)

# ! temp
idx = (500, 833)
b_ = b_field[idx[0], idx[1]] * 1e3
print(f"b at {idx}: {b_}")

idx = (833, 500)
b_ = b_field[idx[0], idx[1]] * 1e3
print(f"b at {idx}: {b_}")

path_r = outer_diameter_m / 2 + 2.5e-3

plot_b_field(b_field, x_grid, y_grid, path_r)

bs = circle_bs(b_field, path_r, grid_n)

max_b_tan_mT = max_b_tan_T(bs) * 1000
max_b_rad_mT = max_b_rad_T(bs) * 1000

print(f"max_x_b_mT: {max_b_tan_mT:.2f} mT")
print(f"max_y_b_mT: {max_b_rad_mT:.2f} mT")
print(f"min_xy_b_mT: {min_b_xy_T(bs) * 1000:.2f} mT")

k = max_b_rad_mT / max_b_tan_mT
print(f"k: {k:.2f}")

plot_path_b_xy(bs)

b at (500, 833): [511.19691755   0.90795166   0.        ]
b at (833, 500): [-94.40547736   0.90795166   0.        ]

max_x_b_mT: 1253.07 mT
max_y_b_mT: -354.79 mT
min_xy_b_mT: 94.41 mT
k: -0.28
bxs: [-94.40547736230057, -94.35648186773163, -94.2422297265605, -94.06289046917678, -93.81873021334489, -93.51011128038337, -93.13749167001784, -92.70142439230268, 1253.0671997449354, -93.40820313035388, -92.77042958252558, 

[cwahn]

I am sorry. The problem was related with the indexing of the grid, not Magpylib. Closing the issue.

[OrtnerMichael]

hey @cwahn

sorry for this slow reaction on our side

Very nice plots you are making :)

one comment: when you have a problem its best to publish only a minimal example that reproduces the error. You posted more than 200 lines of code - its quite readable but still :).

glad you figured it out right away by yourself

br michael