Сегментированный цилиндр Хальбаха: поле и дискретизация намагничивания

"""Самостоятельное полевое исследование во встроенном Python FEMM 5.

Все размеры геометрии указаны в миллиметрах, результаты — в СИ. Значения
гармонических токов являются амплитудами. Материалы и условия имеют русские
имена; номера групп описаны в README. Программа создаёт новую calculated-NN
рядом с собой, не заменяя предшествующих расчётов. Каждая входная модель
сохраняется, закрывается, открывается повторно и только затем решается.
Открытые пользовательские документы не изменяются. Параметры — в CONFIG.
"""

import os
import sys
import json
import math

MU0 = 4e-7 * math.pi
script = globals().get("__file__", "") or (sys.argv[0] if sys.argv else "")
HERE = os.path.dirname(os.path.abspath(script)) if script else os.getcwd()
number = 1
while True:
    OUT = HERE + "/calculated-%02d" % number
    try:
        os.stat(OUT)
        number += 1
    except OSError:
        os.mkdir(OUT)
        break


def scalar(value):
    if isinstance(value, (list, tuple)):
        return scalar(value[0])
    return complex(value).real


def complex_parts(value):
    value = complex(value)
    return [value.real, value.imag]


def point(radius, angle):
    angle = angle * math.pi / 180
    return (radius * math.cos(angle), radius * math.sin(angle))


def line(a, b, group=0, size=0):
    mi_addnode(a[0], a[1])
    mi_addnode(b[0], b[1])
    mi_addsegment(a[0], a[1], b[0], b[1])
    mi_selectsegment((a[0] + b[0]) / 2, (a[1] + b[1]) / 2)
    boundary = (
        "Внешняя граница A=0"
        if globals().get("AXIS", False) and abs(a[0]) + abs(b[0]) < 1e-9
        else "<None>"
    )
    mi_setsegmentprop(boundary, size, 0 if size else 1, 0, group)
    mi_clearselected()


def polygon(points, group=0, size=0):
    for index in range(len(points)):
        line(points[index], points[(index + 1) % len(points)], group, size)


def rectangle(x1, y1, x2, y2, group=0, size=0):
    polygon([(x1, y1), (x2, y1), (x2, y2), (x1, y2)], group, size)


def arc(radius, start, end, group=0, step=2):
    a, b = point(radius, start), point(radius, end)
    mi_addnode(a[0], a[1])
    mi_addnode(b[0], b[1])
    mi_addarc(a[0], a[1], b[0], b[1], end - start, step)
    middle = point(radius, (start + end) / 2)
    mi_selectarcsegment(middle[0], middle[1])
    mi_setarcsegmentprop(step, "<None>", 0, group)
    mi_clearselected()


def circle(radius, group=0, step=2):
    for start in (0, 90, 180, 270):
        arc(radius, start, start + 90, group, step)


def label(x, y, material, mesh, group=0, circuit="<None>", turns=1, angle=0):
    mi_addblocklabel(x, y)
    mi_selectlabel(x, y)
    mi_setblockprop(material, 0, mesh, circuit, angle, group, turns)
    mi_clearselected()


def begin(mesh, kind="planar", frequency=0, depth=40, radius=100, br=1.1, mur=1.05):
    global AXIS
    AXIS = kind == "axi"
    newdocument(0)
    mi_probdef(frequency, "millimeters", kind, 1e-9, depth, 25)
    mi_addmaterial("Воздух", 1, 1, 0, 0, 0)
    mi_addmaterial("Обмотка без вытеснения тока", 1, 1, 0, 0, 0)
    mi_addmaterial("Медь 58 МСм на метр", 1, 1, 0, 0, 58)
    mi_addmaterial("Сталь с линейной проницаемостью", 1000, 1000, 0, 0, 0)
    mi_addmaterial("Сталь с учебной кривой насыщения", 1, 1, 0, 0, 0)
    for b, h in (
        (0, 0),
        (0.3, 100),
        (0.6, 200),
        (1, 350),
        (1.3, 700),
        (1.5, 1800),
        (1.65, 5000),
        (1.8, 15000),
        (2, 60000),
        (2.2, 200000),
    ):
        mi_addbhpoint("Сталь с учебной кривой насыщения", b, h)
    mi_addmaterial("Постоянный магнит учебный NdFeB", mur, mur, br / (MU0 * mur), 0, 0)
    mi_addboundprop("Внешняя граница A=0", 0, 0, 0, 0, 0, 0, 0, 0, 0)
    if kind == "axi":
        rectangle(0, -radius, radius, radius)
        for x, y in ((0, 0), (radius, 0), (radius / 2, -radius), (radius / 2, radius)):
            mi_selectsegment(x, y)
            mi_setsegmentprop("Внешняя граница A=0", 0, 1, 0, 0)
            mi_clearselected()
    else:
        circle(radius, step=4)
        for angle in (45, 135, 225, 315):
            x, y = point(radius, angle)
            mi_selectarcsegment(x, y)
            mi_setarcsegmentprop(4, "Внешняя граница A=0", 0, 0)
            mi_clearselected()


def band(r1, r2, count, fraction, group, mesh, phase=0, mode="radial"):
    """Кольцевая сборка: заданные магнитные секторы и воздушные промежутки."""
    pitch = 360 / count
    for index in range(count):
        start = index * pitch + phase
        stop = start + pitch * fraction
        for radius in (r1, r2):
            arc(radius, start, stop, group, step=2)
            if fraction < 0.999999:
                arc(radius, stop, start + pitch, group, step=2)
        line(point(r1, start), point(r2, start), group)
        line(point(r1, stop), point(r2, stop), group)
        angle = (start + stop) / 2
        x, y = point((r1 + r2) / 2, angle)
        if mode == "steel":
            label(x, y, "Сталь с линейной проницаемостью", mesh, group)
        else:
            direction = 2 * angle if mode == "halbach" else angle + 180 * (index % 2)
            label(x, y, "Постоянный магнит учебный NdFeB", mesh, group, angle=direction)
        if fraction < 0.999999:
            x, y = point((r1 + r2) / 2, (stop + start + pitch) / 2)
            label(x, y, "Воздух", mesh, group)


def energy():
    mo_clearblock()
    mo_groupselectblock()
    values = [scalar(mo_blockintegral(k)) for k in (2, 17, 4)]
    mo_clearblock()
    return values


def solve(tag, group=2, probe=(1, 0), circuits=()):
    path = OUT + "/" + tag + ".fem"
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    e, w, p = energy()
    row = {
        "state": tag,
        "elements": mo_numelements(),
        "energy_J": e,
        "coenergy_J": w,
        "loss_W": p,
    }
    if group:
        mo_groupselectblock(group)
        if CONFIG.get("kind") != "course":
            row["Fx_N"] = scalar(mo_blockintegral(18))
            row["Fy_N"] = scalar(mo_blockintegral(19))
        row["torque_Nm"] = scalar(mo_blockintegral(22))
        if CONFIG.get("kind") != "course":
            row["lorentz_x_N"] = scalar(mo_blockintegral(11))
            row["lorentz_y_N"] = scalar(mo_blockintegral(12))
        mo_clearblock()
    if CONFIG.get("kind") == "gear":
        for part in (2, 3, 4):
            mo_groupselectblock(part)
            row["T%d_Nm" % part] = scalar(mo_blockintegral(22))
            mo_clearblock()
    values = mo_getpointvalues(probe[0], probe[1])
    row["Bx_T"] = scalar(values[1])
    row["By_T"] = scalar(values[2])
    row["Hx_Am"] = scalar(values[5])
    row["Hy_Am"] = scalar(values[6])
    for circuit in circuits:
        i, v, psi = mo_getcircuitproperties(circuit)
        row[circuit + "_I_A"] = scalar(i)
        row[circuit + "_Psi_Wb"] = scalar(psi)
        row[circuit + "_Psi_im_Wb"] = complex(psi).imag
        row[circuit + "_V_V"] = scalar(v)
        row[circuit + "_V_im_V"] = complex(v).imag
    mo_close()
    mi_close()
    return row


def check(name, value, limit, passed):
    CHECKS.append(
        {"name": name, "value": value, "limit": limit, "passed": bool(passed)}
    )


def relative(a, b, minimum=1e-15):
    return abs(a - b) / max(abs(a), abs(b), minimum)


def finish(plot=None):
    result = {
        "title": CONFIG["title"],
        "records": ROWS,
        "checks": CHECKS,
        "units": "Геометрия мм; магнитное поле Тл, А/м; энергия Дж; сила Н; момент Н·м; потокосцепление Вб.",
        "passed": bool(CHECKS) and all(row["passed"] for row in CHECKS),
    }
    if plot:
        result["plot"] = plot
    with open(OUT + "/result.json", "w") as stream:
        json.dump(result, stream)
    print(CONFIG["id"], result["passed"], len(ROWS), "решений", OUT)
    assert result["passed"], "Не пройдена независимая проверка физической модели."


ROWS = []
CHECKS = []


CONFIG = {
    "id": "halbach-cylinder",
    "kind": "halbach",
    "title": "Сегментированный цилиндр Хальбаха: поле и дискретизация намагничивания",
    "baseline": 16,
    "metric": "Bx_T",
    "sweep": [8, 16, 32],
    "parameter_name": "Число магнитных сегментов",
    "metric_label": "Центральная индукция, Тл",
    "dimensions": "Внутренний радиус 15 мм, наружный 30 мм, Br=1,1 Тл, μr=1, глубина 40 мм, внешняя граница радиуса 150 мм.",
}


def build(mesh, parameter, current=1):
    kind = CONFIG["kind"]
    axis = kind in (
        "plunger",
        "coil-magnet",
        "pm-bearing",
        "sphere-temperature",
        "clutch",
        "levitation",
    )
    frequency = parameter if kind == "levitation" else 0
    br = 1.1 * (1 - 0.001 * (parameter - 20)) if kind == "sphere-temperature" else 1.1
    begin(
        mesh,
        "axi" if axis else "planar",
        frequency,
        40,
        150 if kind == "halbach" else 100,
        br,
        1 if kind == "halbach" else 1.05,
    )
    coil = "Обмотка без вытеснения тока"
    steel = "Сталь с линейной проницаемостью"
    magnet = "Постоянный магнит учебный NdFeB"
    circuits = []
    if kind == "halbach":
        band(15, 30, int(parameter), 1, 2, mesh, mode="halbach")
        label(0, 0, "Воздух", mesh)
        label(70, 0, "Воздух", 6 * mesh)
        probe = (0, 0)
    elif kind == "gear":
        circle(10, 2)
        label(0, 0, "Воздух", 4 * mesh, 2)
        label(13, 0, steel, mesh, 2)
        band(16, 22, 4, 0.94, 2, mesh, parameter)
        band(24, 28, 8, 0.5, 4, mesh, mode="steel")
        band(30, 36, 12, 0.94, 3, mesh)
        circle(42, 3)
        label(39, 0, steel, mesh, 3)
        label(23, 0, "Воздух", mesh)
        label(29, 0, "Воздух", mesh)
        label(70, 0, "Воздух", 6 * mesh)
        probe = (23, 0)
    elif kind == "reluctance-rotor":
        points = []
        for index in range(360):
            theta = index * 2 * math.pi / 360
            x, y = 35 * math.cos(theta), 20 * math.sin(theta)
            turn = parameter * math.pi / 180
            points.append(
                (
                    x * math.cos(turn) - y * math.sin(turn),
                    x * math.sin(turn) + y * math.cos(turn),
                )
            )
        polygon(points, 2)
        label(0, 0, steel, mesh, 2)
        circle(42)
        label(40, 0, "Воздух", 0.6 * mesh)
        mi_addcircprop("Эквивалентная фазная обмотка", current, 1)
        for x1, x2, sign in ((-60, -50, 1), (50, 60, -1)):
            rectangle(x1, -25, x2, 25)
            label(
                (x1 + x2) / 2,
                0,
                coil,
                1.5 * mesh,
                1,
                "Эквивалентная фазная обмотка",
                200 * sign,
            )
        label(80, 0, "Воздух", 5 * mesh)
        circuits = ["Эквивалентная фазная обмотка"]
        probe = (0, 0)
    elif kind in ("reactor", "nonlinear-reactor", "transformer-matrix"):
        gap = parameter if kind == "reactor" else 0.5
        material = (
            "Сталь с учебной кривой насыщения" if kind == "nonlinear-reactor" else steel
        )
        polygon(
            [
                (-35, -30),
                (35, -30),
                (35, -gap / 2),
                (25, -gap / 2),
                (25, -20),
                (-25, -20),
                (-25, 20),
                (25, 20),
                (25, gap / 2),
                (35, gap / 2),
                (35, 30),
                (-35, 30),
            ],
            size=min(0.25, gap / 3) * mesh,
        )
        label(-30, 0, material, 2 * mesh, 1)
        label(0, 0, "Воздух", 5 * mesh)
        mi_addcircprop("Первичная обмотка", current, 1)
        for x1, x2, sign in ((-43, -37, 1), (-23, -17, -1)):
            rectangle(x1, -14, x2, 14)
            label(
                (x1 + x2) / 2, 0, coil, 1.5 * mesh, 0, "Первичная обмотка", 100 * sign
            )
        circuits = ["Первичная обмотка"]
        if kind == "transformer-matrix":
            mi_addcircprop("Вторичная обмотка", parameter, 1)
            for x1, x2, sign in ((17, 23, -1), (37, 43, 1)):
                rectangle(x1, -14, x2, 14)
                label(
                    (x1 + x2) / 2,
                    0,
                    coil,
                    1.5 * mesh,
                    0,
                    "Вторичная обмотка",
                    50 * sign,
                )
            circuits.append("Вторичная обмотка")
        probe = (30, 0)
    elif kind == "plunger":
        mi_addcircprop("Катушка привода", current, 1)
        rectangle(10, -10, 16, 10)
        label(13, 0, coil, mesh, 1, "Катушка привода", 200)
        rectangle(0, parameter - 8, 8, parameter + 8, 2)
        label(4, parameter, steel, mesh, 2)
        label(35, 0, "Воздух", 5 * mesh)
        circuits = ["Катушка привода"]
        probe = (4, parameter)
    elif kind == "coil-magnet":
        mi_addcircprop("Катушка привода", current, 1)
        rectangle(10, -6, 14, 6)
        label(12, 0, coil, mesh, 1, "Катушка привода", 100)
        rectangle(0, parameter - 3, 6, parameter + 3, 2)
        label(3, parameter, magnet, 0.2 * mesh, 2, angle=90)
        rectangle(0, parameter - 7, 9, parameter + 7)
        label(8, parameter, "Воздух", 0.25 * mesh)
        label(35, 0, "Воздух", 5 * mesh)
        circuits = ["Катушка привода"]
        probe = (3, parameter)
    elif kind == "pm-bearing":
        for center, angle, group in ((0, 90, 1), (parameter + 4, 270, 2)):
            rectangle(5, center - 2, 12, center + 2, group)
            label(8, center, magnet, 0.3 * mesh, group, angle=angle)
        rectangle(0, -5, 16, parameter + 9)
        label(14, 0, "Воздух", 0.35 * mesh)
        label(40, 0, "Воздух", 5 * mesh)
        probe = (8, 2 + parameter / 2)
    elif kind == "sphere-temperature":
        mi_addnode(0, -5)
        mi_addnode(0, 5)
        mi_addarc(0, -5, 0, 5, 180, max(0.25, mesh))
        line((0, -5), (0, 5), 2)
        label(2, 0, magnet, mesh / 5, 2, angle=90)
        mi_addnode(0, -15)
        mi_addnode(0, 15)
        mi_addarc(0, -15, 0, 15, 180, max(0.25, mesh))
        label(10, 0, "Воздух", 0.4 * mesh)
        label(30, 0, "Воздух", 5 * mesh)
        probe = (1, 0)
    elif kind == "clutch":
        polygon(
            [
                (0, -15),
                (30, -15),
                (30, 0),
                (25, 0),
                (25, -10),
                (10, -10),
                (10, 0),
                (0, 0),
            ],
            1,
            size=mesh,
        )
        label(5, -5, steel, 2 * mesh, 1)
        rectangle(11, -9, 24, -1)
        mi_addcircprop("Катушка муфты", current, 1)
        label(17, -5, coil, mesh, 1, "Катушка муфты", 200)
        rectangle(
            0, parameter, 30, parameter + 5, 2, size=min(0.3, parameter / 3) * mesh
        )
        label(15, parameter + 2, steel, mesh, 2)
        label(40, 5, "Воздух", 5 * mesh)
        circuits = ["Катушка муфты"]
        probe = (5, parameter / 2)
    elif kind == "levitation":
        mi_addcircprop("Катушка переменного тока", current, 1)
        rectangle(8, -10, 20, -2)
        label(14, -6, coil, mesh, 1, "Катушка переменного тока", 100)
        rectangle(0, 3, 25, 6, 2, size=mesh)
        label(12, 4.5, "Медь 58 МСм на метр", mesh / 2, 2)
        label(45, 0, "Воздух", 5 * mesh)
        circuits = ["Катушка переменного тока"]
        probe = (12, 4.5)
    elif kind == "linear-motor":
        for index in range(8):
            x = -40 + index * 10
            rectangle(x, -10, x + 9, -4, 1)
            label(x + 4, -7, magnet, mesh, 1, angle=90 + 180 * (index % 2))
        mi_addcircprop("Подвижная катушка", current, 1)
        for x1, x2, sign in (
            (parameter - 11, parameter - 6, 1),
            (parameter - 1, parameter + 4, -1),
        ):
            rectangle(x1, 1, x2, 5, 2)
            label((x1 + x2) / 2, 3, coil, 0.2 * mesh, 2, "Подвижная катушка", 50 * sign)
        rectangle(parameter - 13, 0, parameter + 6, 8)
        label(parameter - 12, 0.5, "Воздух", 0.25 * mesh)
        label(60, 20, "Воздух", 5 * mesh)
        circuits = ["Подвижная катушка"]
        probe = (parameter, 0)
    return probe, circuits


kind = CONFIG["kind"]
baseline = CONFIG["baseline"]
metric = CONFIG["metric"]
mesh_rows = []
meshes = CONFIG.get("meshes", (1.2, 0.8, 0.5))
fine_h = meshes[-1]
for index, mesh in enumerate(meshes):
    parameter = 0 if kind == "transformer-matrix" else baseline
    probe, circuits = build(mesh, parameter, CONFIG.get("current", 1))
    row = solve(
        "сетка-%d" % (index + 1),
        0
        if kind
        in (
            "reactor",
            "nonlinear-reactor",
            "transformer-matrix",
            "halbach",
            "sphere-temperature",
        )
        else 2,
        probe,
        circuits,
    )
    row["mesh_mm"] = mesh
    row["parameter"] = (
        CONFIG.get("current", 1) if kind == "nonlinear-reactor" else parameter
    )
    mesh_rows.append(row)
    ROWS.append(row)
fine = mesh_rows[-1]
error = relative(mesh_rows[-2][metric], fine[metric])
check(
    "Стабилизация наблюдаемой величины на двух мелких сетках",
    error,
    0.025,
    error < 0.025,
)
check(
    "Действительное уточнение сетки",
    fine["elements"] / mesh_rows[0]["elements"],
    1.2,
    fine["elements"] > 1.2 * mesh_rows[0]["elements"],
)
if kind == "halbach":
    count = int(baseline)
    exact = 1.1 * math.log(2) * math.sin(2 * math.pi / count) / (2 * math.pi / count)
    err = relative(abs(fine["Bx_T"]), exact)
    check(
        "Поле сегментированного цилиндра и аналитическая гармоника",
        err,
        0.025,
        err < 0.025,
    )
elif kind == "sphere-temperature":
    exact = 2 * 1.1 * (1 - 0.001 * (baseline - 20)) / (1.05 + 2)
    err = relative(fine["By_T"], exact)
    check("Внутреннее поле равномерно намагниченного шара", err, 0.005, err < 0.005)
elif kind in (
    "reactor",
    "nonlinear-reactor",
    "transformer-matrix",
    "plunger",
    "clutch",
    "reluctance-rotor",
):
    name = circuits[0]
    expected = 0.5 * CONFIG.get("current", 1) * fine[name + "_Psi_Wb"]
    if kind != "nonlinear-reactor":
        err = relative(expected, fine["energy_J"])
        check("Энергия поля и работа токового источника", err, 0.002, err < 0.002)
    else:
        check(
            "Положительная энергия при насыщении",
            fine["energy_J"],
            0,
            fine["energy_J"] > 0,
        )

for value in CONFIG.get("sweep", []):
    current = CONFIG.get("current", 1)
    parameter = value
    if kind == "nonlinear-reactor":
        current = value
        parameter = 0.5
    probe, circuits = build(fine_h, parameter, current)
    row = solve(
        "параметр-%02d" % len(ROWS),
        0
        if kind
        in (
            "reactor",
            "nonlinear-reactor",
            "transformer-matrix",
            "halbach",
            "sphere-temperature",
        )
        else 2,
        probe,
        circuits,
    )
    row["parameter"] = value
    row["mesh_mm"] = fine_h
    if kind == "nonlinear-reactor":
        row["static_inductance_H"] = row["Первичная обмотка_Psi_Wb"] / current
    ROWS.append(row)

if kind == "reactor":
    rows = ROWS[3:]
    check(
        "Индуктивность уменьшается при увеличении зазора",
        rows[-1]["energy_J"] / rows[0]["energy_J"],
        1,
        rows[-1]["energy_J"] < rows[0]["energy_J"],
    )
elif kind == "nonlinear-reactor":
    rows = ROWS[3:]
    check(
        "Насыщение снижает статическую индуктивность",
        rows[-1]["static_inductance_H"] / rows[0]["static_inductance_H"],
        1,
        rows[-1]["static_inductance_H"] < rows[0]["static_inductance_H"],
    )
elif kind == "transformer-matrix":
    probe, circuits = build(fine_h, 1, 0)
    secondary = solve("вторичное-возбуждение", 0, probe, circuits)
    ROWS.append(secondary)
    m12 = secondary["Первичная обмотка_Psi_Wb"]
    m21 = fine["Вторичная обмотка_Psi_Wb"] / CONFIG.get("current", 1)
    err = relative(m12, m21)
    check("Взаимность индуктивностей", err, 0.005, err < 0.005)
    l1 = fine["Первичная обмотка_Psi_Wb"] / CONFIG.get("current", 1)
    l2 = secondary["Вторичная обмотка_Psi_Wb"]
    k = abs(m12) / math.sqrt(l1 * l2)
    check("Пассивность матрицы индуктивности", k, 1, k < 1 and k > 0)
elif kind in ("plunger", "clutch"):
    controls = []
    for value in (baseline - 0.2, baseline + 0.2):
        probe, circuits = build(fine_h, value, CONFIG.get("current", 1))
        r = solve("виртуальное-смещение-%d" % len(ROWS), 2, probe, circuits)
        ROWS.append(r)
        controls.append(r)
    derivative = (controls[1]["coenergy_J"] - controls[0]["coenergy_J"]) / 0.0004
    err = relative(derivative, fine["Fy_N"])
    check("Сила по коэнергии и взвешенному тензору напряжений", err, 0.06, err < 0.06)
elif kind == "reluctance-rotor":
    controls = []
    for value in (baseline - 0.5, baseline + 0.5):
        probe, circuits = build(fine_h, value, CONFIG.get("current", 1))
        r = solve("виртуальный-поворот-%d" % len(ROWS), 2, probe, circuits)
        ROWS.append(r)
        controls.append(r)
    derivative = (controls[1]["coenergy_J"] - controls[0]["coenergy_J"]) / (
        math.pi / 180
    )
    err = relative(derivative, fine["torque_Nm"])
    check("Момент и угловая производная коэнергии", err, 0.03, err < 0.03)
elif kind in ("coil-magnet", "linear-motor"):
    probe, circuits = build(fine_h, baseline, -CONFIG.get("current", 1))
    reversed_row = solve("обратный-ток", 2, probe, circuits)
    ROWS.append(reversed_row)
    field = "Fx_N" if kind == "linear-motor" else "Fy_N"
    err = relative(fine[field], -reversed_row[field])
    check("Смена направления силы при обращении тока", err, 0.01, err < 0.01)
    if kind == "linear-motor":
        err = relative(fine["Fx_N"], fine["lorentz_x_N"])
        check("Сила Лоренца и взвешенный тензор напряжений", err, 0.04, err < 0.04)
elif kind == "pm-bearing":
    rows = ROWS[3:]
    check(
        "Репульсивная осевая сила",
        min(r["Fy_N"] for r in rows),
        0,
        all(r["Fy_N"] > 0 for r in rows),
    )
    check(
        "Спад силы при увеличении зазора",
        rows[-1]["Fy_N"] / rows[0]["Fy_N"],
        1,
        rows[-1]["Fy_N"] < rows[0]["Fy_N"],
    )
elif kind == "levitation":
    rows = ROWS[3:]
    check(
        "Отсутствие вихревых потерь при постоянном токе",
        rows[0]["loss_W"],
        1e-8,
        abs(rows[0]["loss_W"]) < 1e-8,
    )
    check(
        "Положительные вихревые потери",
        min(r["loss_W"] for r in rows[1:]),
        0,
        all(r["loss_W"] > 0 for r in rows[1:]),
    )
    check("Сила отталкивания проводящего диска", fine["Fy_N"], 0, fine["Fy_N"] > 0)
elif kind == "gear":
    balance = abs(fine["T2_Nm"] + fine["T3_Nm"] + fine["T4_Nm"]) / max(
        abs(fine["T2_Nm"]), abs(fine["T3_Nm"])
    )
    check("Баланс моментов трёх магнитных тел", balance, 0.05, balance < 0.05)
    a = 0j
    b = 0j
    for row in ROWS[3:]:
        phase = 2 * row["parameter"] * math.pi / 180
        factor = complex(math.cos(phase), -math.sin(phase))
        a += row["T2_Nm"] * factor
        b += row["T3_Nm"] * factor
    ratio = abs(a / b)
    check(
        "Отношение основных гармоник моментов после отделения зубцовых составляющих",
        ratio,
        1 / 3,
        relative(ratio, 1 / 3) < 0.05,
    )

finish(
    {
        "x": "parameter",
        "y": metric,
        "x_label": CONFIG.get("parameter_name", "Параметр"),
        "y_label": CONFIG.get("metric_label", metric),
    }
)