ЭДС машины с постоянными магнитами: угол, потокосцепление и скорость

"""Самостоятельное полевое исследование во встроенном 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": "pm-machine-emf",
    "kind": "course",
    "study": "pm-emf",
    "title": "ЭДС машины с постоянными магнитами: угол, потокосцепление и скорость",
    "source_model": "spm",
    "baseline": 15,
    "angle": 15,
    "metric": "Фаза A_Psi_Wb",
    "sweep": [
        0.0,
        7.5,
        15.0,
        22.5,
        30.0,
        37.5,
        45.0,
        52.5,
        60.0,
        67.5,
        75.0,
        82.5,
        90.0,
        97.5,
        105.0,
        112.5,
        120.0,
        127.5,
        135.0,
        142.5,
        150.0,
        157.5,
        165.0,
        172.5,
        180.0,
    ],
    "metric_label": "Потокосцепление фазы A, Вб",
    "parameter_name": "Механический угол, град",
    "circuits": ["Фаза A", "Фаза B", "Фаза C"],
    "group": 0,
    "probe": [58.5, 0],
    "intro": "Магниты ротора создают периодическое потокосцепление неподвижной фазной обмотки. Из последовательности статических полей получают ЭДС при заданной скорости; сама магнитостатическая задача не содержит вращения во времени.",
    "theory": "Закон Фарадея связывает ЭДС с производной полного потокосцепления. Здесь механический угол задаётся в радианах при дифференцировании и в градусах только при построении геометрии. Число пар полюсов равно двум, поэтому половина механического оборота содержит один электрический период. Периодический ряд не должен дважды включать крайний отсчёт в центральной разности. Масштабирование скорости меняет ЭДС, но не исходную магнитостатическую картину при нулевом токе. Используется периодическая центральная разность четвёртого порядка и шаг 7,5° механического угла. Основная гармоника ЭДС дополнительно сравнивается с обработкой каждого второго полевого положения, то есть с шагом 15°. Эта угловая проверка отделена от проверки пространственной сетки; точность высоких гармоник требует своего более мелкого углового ряда.",
    "formula": "e_A=-\\omega_m\\,d\\Psi_A/d\\theta_m,\\qquad f=pn/60",
    "scope": "Двумерная учебная геометрия; точные входы и русские свойства включены в комплект. Торцевые поля, механическая динамика, температурная зависимость и реальные потери стали требуют отдельных моделей.",
    "related_lessons": ["I10.1", "II07.2", "meas-search-coil"],
    "dimensions": "Отдельная русская учебная копия геометрии нашего курса. Все координаты, витки, материалы и глубина заданы в input.fem; рабочие токи и угловая серия задаются программой. Сталь линейная, μr=1000.",
    "plot_metric": "emf_A_V",
    "plot_label": "ЭДС фазы A при 1000 об/мин, В",
    "mesh_scales": [1, 0.7, 0.5],
    "magnet_hc": 909456.8176679733,
    "depth": 80.0,
}

with open(HERE + "/input.fem", "r") as stream:
    BASE_TEXT = stream.read()


def prepare(scale=1, angle=0, currents=None, frequency=0, defect=False):
    lines = BASE_TEXT.split("\n")
    start = next(
        i for i in range(len(lines)) if lines[i].startswith("[NumBlockLabels]")
    )
    count = int(lines[start].split("=")[1])
    for i in range(start + 1, start + 1 + count):
        fields = lines[i].split("\t")
        fields[3] = str(float(fields[3]) * scale)
        lines[i] = "\t".join(fields)
    path = OUT + "/текущая-настройка.fem"
    with open(path, "w") as stream:
        stream.write("\n".join(lines) + "\n")
    opendocument(path)
    mi_probdef(frequency, "millimeters", "planar", 1e-9, CONFIG.get("depth", 80), 25)
    for name, value in (currents or {}).items():
        mi_modifycircprop(name, 1, value)
    if defect:
        mi_addmaterial(
            "NdFeB с ослабленной намагниченностью",
            1.05,
            1.05,
            CONFIG["magnet_hc"] * 0.6,
            0,
            0,
        )
        mi_selectlabel(52, 0)
        mi_setblockprop(
            "NdFeB с ослабленной намагниченностью",
            0,
            1.5 * scale,
            "<None>",
            "theta",
            2,
            1,
        )
        mi_clearselected()
    if angle:
        mi_selectgroup(2)
        mi_moverotate(0, 0, angle)
        mi_clearselected()


def phase_currents(angle=0, amplitude=1, delta=0):
    return {
        name: amplitude * math.cos((2 * angle + delta - offset) * math.pi / 180)
        for name, offset in (("Фаза A", 0), ("Фаза B", 120), ("Фаза C", 240))
    }


study = CONFIG["study"]
names = CONFIG["circuits"]
opendocument(HERE + "/input.fem")
ROWS.append(solve("input", CONFIG.get("group", 0), CONFIG.get("probe", [1, 0]), names))
mesh_rows = []
mesh_scales = CONFIG.get("mesh_scales", (0.8, 0.5, 0.3))
fine_scale = mesh_scales[-1]
for index, scale in enumerate(mesh_scales):
    angle = CONFIG.get("angle", 0)
    if study in ("pm-emf", "pm-defect"):
        currents = {name: 0 for name in names}
    elif study == "pm-skew":
        currents = phase_currents(angle, 3, 60)
    elif study == "dc-reaction":
        currents = {names[0]: 1, names[1]: 2}
    elif study == "ct-burden":
        currents = {names[0]: 1, names[1]: 0}
    elif study == "three-phase":
        currents = {names[0]: 1, names[1]: 0, names[2]: 0}
    elif study == "induction-slip":
        currents = {
            names[0]: 1,
            names[1]: complex(-0.5, -math.sqrt(3) / 2),
            names[2]: complex(-0.5, math.sqrt(3) / 2),
        }
    else:
        currents = {
            names[0]: 1,
            names[1]: -0.5 if study == "synrm-angle" else 0,
            names[2]: -0.5 if study == "synrm-angle" else 0,
        }
    prepare(
        scale,
        angle,
        currents,
        frequency=CONFIG.get("frequency", 0),
        defect=study == "pm-defect",
    )
    row = solve(
        "сетка-%d" % (index + 1),
        CONFIG.get("group", 0),
        CONFIG.get("probe", [1, 0]),
        names,
    )
    row["mesh_scale"] = scale
    row["parameter"] = (
        2
        if study == "dc-reaction"
        else CONFIG.get("frequency", 0)
        if study == "induction-slip"
        else angle
    )
    ROWS.append(row)
    mesh_rows.append(row)
fine = mesh_rows[-1]
metric = CONFIG["metric"]
err = relative(mesh_rows[-2][metric], fine[metric])
check("Стабилизация наблюдаемой величины при сгущении сетки", err, 0.025, err < 0.025)
check(
    "Фактическое увеличение числа элементов",
    fine["elements"] / mesh_rows[0]["elements"],
    1.2,
    fine["elements"] > 1.2 * mesh_rows[0]["elements"],
)

if study in ("pm-emf", "pm-skew", "pm-defect", "synrm-angle", "srm-map"):
    for angle in CONFIG["sweep"]:
        if study in ("pm-emf", "pm-defect"):
            currents = {name: 0 for name in names}
        elif study == "pm-skew":
            currents = phase_currents(angle, 3, 60)
        else:
            currents = {
                names[0]: 1,
                names[1]: -0.5 if study == "synrm-angle" else 0,
                names[2]: -0.5 if study == "synrm-angle" else 0,
            }
        prepare(fine_scale, angle, currents, defect=study == "pm-defect")
        row = solve(
            "угол-" + str(angle),
            CONFIG.get("group", 0),
            CONFIG.get("probe", [1, 0]),
            names,
        )
        row["parameter"] = angle
        row["mesh_scale"] = fine_scale
        ROWS.append(row)
    data = ROWS[4:]
    if study in ("pm-emf", "pm-skew", "pm-defect"):
        # Первый и последний узлы — одинаковое положение после периода.
        key = names[0] + "_Psi_Wb" if study != "pm-skew" else "torque_Nm"
        amplitude = max(abs(r[key]) for r in data)
        err = abs(data[-1][key] - data[0][key]) / max(amplitude, 1e-12)
        check("Замыкание углового периода", err, 0.025, err < 0.025)
        if study in ("pm-emf", "pm-defect"):
            omega = 1000 * 2 * math.pi / 60
            step = (data[1]["parameter"] - data[0]["parameter"]) * math.pi / 180
            values = [r[names[0] + "_Psi_Wb"] for r in data[:-1]]

            def derivative_emf(values, step):
                size = len(values)
                return [
                    -omega
                    * (
                        -values[(i + 2) % size]
                        + 8 * values[(i + 1) % size]
                        - 8 * values[(i - 1) % size]
                        + values[(i - 2) % size]
                    )
                    / (12 * step)
                    for i in range(size)
                ]

            emf = derivative_emf(values, step)
            coarse_emf = derivative_emf(values[::2], 2 * step)
            order = int(round((CONFIG["sweep"][-1] - CONFIG["sweep"][0]) * 2 / 360))

            def amplitude(samples):
                size = len(samples)
                coefficient = 0j
                for i, value in enumerate(samples):
                    phase = 2 * math.pi * order * i / size
                    coefficient += value * complex(math.cos(phase), -math.sin(phase))
                return 2 * abs(coefficient) / size

            fine_amplitude = amplitude(emf)
            coarse_amplitude = amplitude(coarse_emf)
            angular_error = relative(fine_amplitude, coarse_amplitude)
            check(
                "Сходимость основной гармоники ЭДС по угловому шагу 15 и 7,5 град",
                angular_error,
                0.01,
                angular_error < 0.01,
            )
            for i, value in enumerate(emf):
                data[i]["emf_A_V"] = value
            check(
                "Ненулевая ЭДС при вращении",
                max(abs(v) for v in emf),
                0,
                max(abs(v) for v in emf) > 0,
            )
        else:
            values = []
            averaged = []
            for center in data:
                angle = center["parameter"]
                currents = phase_currents(angle, 3, 60)
                torques = [center["torque_Nm"]]
                for offset in (-5, 5):
                    # У всех сечений одного положения ОДИН И ТОТ ЖЕ ток статора.
                    prepare(fine_scale, angle + offset, currents)
                    row = solve(
                        "сечение-%03d-%s" % (angle, str(offset)),
                        2,
                        CONFIG.get("probe", [1, 0]),
                        names,
                    )
                    row.update(global_angle_deg=angle, slice_offset_deg=offset)
                    ROWS.append(row)
                    torques.append(row["torque_Nm"])
                center["skew_averaged_torque_Nm"] = sum(torques) / 3
                if angle != CONFIG["sweep"][-1]:
                    values.append(center["torque_Nm"])
                    averaged.append(sum(torques) / 3)
            factor = (1 + 2 * math.cos(2 * 5 * math.pi / 180)) / 3
            err = relative(sum(averaged), factor * sum(values))
            check(
                "Средний момент и гармонический коэффициент трёх фиксированных токовых сечений",
                err,
                0.03,
                err < 0.03,
            )
    elif study == "srm-map":
        angle = CONFIG.get("angle", 15)
        prepare(fine_scale, angle, {names[0]: 2, names[1]: 0, names[2]: 0})
        double = solve("удвоенный-ток", 0, CONFIG.get("probe", [1, 0]), names)
        ROWS.append(double)
        err = relative(double[names[0] + "_Psi_Wb"], 2 * fine[names[0] + "_Psi_Wb"])
        check("Линейное масштабирование потокосцепления", err, 0.001, err < 0.001)
        err = relative(double["coenergy_J"], 4 * fine["coenergy_J"])
        check("Квадратичное масштабирование коэнергии", err, 0.001, err < 0.001)
    else:
        controls = []
        for angle in (CONFIG["angle"] - 0.5, CONFIG["angle"] + 0.5):
            prepare(0.3, angle, {names[0]: 1, names[1]: -0.5, names[2]: -0.5})
            row = solve(
                "виртуальное-смещение-%d" % len(ROWS),
                2,
                CONFIG.get("probe", [1, 0]),
                names,
            )
            ROWS.append(row)
            controls.append(row)
        derivative = (controls[1]["coenergy_J"] - controls[0]["coenergy_J"]) / (
            math.pi / 180
        )
        err = relative(derivative, fine["torque_Nm"])
        check(
            "Момент и производная коэнергии при фиксированных токах",
            err,
            0.07,
            err < 0.07,
        )
elif study == "dc-reaction":
    separate = []
    for field, armature in ((1, 0), (0, 2)):
        prepare(fine_scale, 0, {names[0]: field, names[1]: armature})
        row = solve("раздельное-возбуждение-%d" % len(ROWS), 0, CONFIG["probe"], names)
        ROWS.append(row)
        separate.append(row)
    err = max(
        relative(fine[name + "_Psi_Wb"], sum(r[name + "_Psi_Wb"] for r in separate))
        for name in names
    )
    check(
        "Суперпозиция потокосцеплений двух систем возбуждения", err, 0.002, err < 0.002
    )
    for armature in (0, 1, 2, 3, 4):
        prepare(fine_scale, 0, {names[0]: 1, names[1]: armature})
        row = solve("реакция-якоря-%d" % armature, 0, CONFIG["probe"], names)
        row["parameter"] = armature
        ROWS.append(row)
elif study == "ct-burden":
    prepare(fine_scale, 0, {names[0]: 0, names[1]: 1})
    secondary = solve("вторичный-столбец", 0, CONFIG["probe"], names)
    ROWS.append(secondary)
    l1 = fine[names[0] + "_Psi_Wb"]
    m1 = fine[names[1] + "_Psi_Wb"]
    m2 = secondary[names[0] + "_Psi_Wb"]
    l2 = secondary[names[1] + "_Psi_Wb"]
    err = relative(m1, m2)
    check("Взаимность связи обмоток", err, 0.001, err < 0.001)
    for resistance in (0.01, 0.03, 0.1, 0.3, 1):
        omega = 2 * math.pi * 50
        ip = 100
        secondary_current = -1j * omega * m1 * ip / (resistance + 1j * omega * l2)
        prepare(fine_scale, 0, {names[0]: ip, names[1]: secondary_current}, 50)
        row = solve("нагрузка-%d" % len(ROWS), 0, CONFIG["probe"], names)
        i2 = secondary_current
        v2 = complex(row[names[1] + "_V_V"], row[names[1] + "_V_im_V"])
        v1 = complex(row[names[0] + "_V_V"], row[names[0] + "_V_im_V"])
        residual = abs(v2 + resistance * i2) / max(abs(v2), 1e-12)
        supplied = 0.5 * (v1 * ip).real
        load = 0.5 * resistance * abs(i2) ** 2
        err = relative(supplied, load)
        row.update(
            parameter=resistance,
            secondary_amplitude_A=abs(i2),
            load_power_W=load,
            supplied_power_W=supplied,
            kvl_residual=residual,
            power_error=err,
        )
        ROWS.append(row)
        check(
            "Невязка уравнения вторичной цепи R=" + str(resistance),
            residual,
            0.002,
            residual < 0.002,
        )
        check("Баланс активной мощности R=" + str(resistance), err, 0.002, err < 0.002)
elif study == "three-phase":
    matrix = []
    for i in range(3):
        prepare(
            fine_scale, 0, {name: 1 if j == i else 0 for j, name in enumerate(names)}
        )
        row = solve("столбец-%d" % i, 0, CONFIG["probe"], names)
        ROWS.append(row)
        row["parameter"] = i
        matrix.append([row[name + "_Psi_Wb"] for name in names])
    asym = max(relative(matrix[i][j], matrix[j][i]) for i in range(3) for j in range(3))
    check("Взаимность трёхфазной матрицы", asym, 0.002, asym < 0.002)
    for label_, currents in (
        ("нулевая-последовательность", [1, 1, 1]),
        ("несимметрия", [1, 2, -1]),
        ("баланс", [1, -0.5, -0.5]),
    ):
        prepare(fine_scale, 0, {name: currents[i] for i, name in enumerate(names)})
        row = solve(label_, 0, CONFIG["probe"], names)
        ROWS.append(row)
        expected = 0.5 * sum(
            currents[i] * matrix[j][i] * currents[j] for i in range(3) for j in range(3)
        )
        err = relative(expected, row["energy_J"])
        check("Энергия поля и матрица: " + label_, err, 0.002, err < 0.002)
elif study == "induction-slip":
    for frequency in CONFIG["sweep"]:
        currents = {
            names[0]: 1,
            names[1]: complex(-0.5, -math.sqrt(3) / 2),
            names[2]: complex(-0.5, math.sqrt(3) / 2),
        }
        prepare(fine_scale, 0, currents, frequency)
        row = solve("частота-%d" % len(ROWS), 0, CONFIG["probe"], names)
        supplied = 0
        for name, current in currents.items():
            voltage = complex(row[name + "_V_V"], row[name + "_V_im_V"])
            supplied += 0.5 * (voltage * complex(current).conjugate()).real
        err = relative(supplied, row["loss_W"])
        row.update(parameter=frequency, terminal_power_W=supplied, power_error=err)
        ROWS.append(row)
        check(
            "Баланс потерь поля и выводов, f=" + str(frequency), err, 0.02, err < 0.02
        )

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