"""Самостоятельное полевое исследование во встроенном 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": "travelling-field-thrust",
    "kind": "linear-induction",
    "title": "Бегущее поле: наведённые токи, тяга пластины и реверс фаз",
}


def travelling(mesh, frequency, sequence=1, amplitude=100, conductivity=35, radius=400):
    begin(mesh, "planar", frequency, 40, radius)
    mi_addmaterial("Алюминиевая вторичная пластина", 1, 1, 0, 0, conductivity)
    mi_addcircprop("Изолированная вторичная пластина", 0, 1)
    rectangle(-180, -20, 180, 40, size=4 * mesh)
    rectangle(-160, -4, 160, 0, 2, size=mesh)
    label(
        0,
        -2,
        "Алюминиевая вторичная пластина",
        mesh,
        2,
        "Изолированная вторичная пластина",
        1,
    )
    circuits = []
    for phase, angle in (("A", 0), ("B", -120), ("C", 120)):
        name = "Фаза " + phase
        radians = sequence * angle * math.pi / 180
        value = amplitude * complex(math.cos(radians), math.sin(radians))
        mi_addcircprop(name, value, 1)
        circuits.append(name)
    for x, phase, sign in (
        (-125, "A", 1),
        (-75, "C", -1),
        (-25, "B", 1),
        (25, "A", -1),
        (75, "C", 1),
        (125, "B", -1),
    ):
        rectangle(x - 10, 10, x + 10, 22, size=mesh)
        label(x, 16, "Обмотка без вытеснения тока", mesh, 1, "Фаза " + phase, 50 * sign)
    label(0, 5, "Воздух", 2 * mesh)
    label(0, 200, "Воздух", 12 * mesh)
    return circuits


def field(
    tag, mesh, frequency=50, sequence=1, amplitude=100, conductivity=35, radius=400
):
    circuits = travelling(mesh, frequency, sequence, amplitude, conductivity, radius)
    path = OUT + "/" + tag + ".fem"
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    mo_groupselectblock(2)
    loss = scalar(mo_blockintegral(4))
    fx = scalar(mo_blockintegral(18))
    fy = scalar(mo_blockintegral(19))
    lorentz = scalar(mo_blockintegral(11))
    plate_current = complex(mo_blockintegral(7))
    mo_clearblock()
    active = 0
    for name in circuits:
        current, voltage, flux = mo_getcircuitproperties(name)
        active += 0.5 * (complex(voltage) * complex(current).conjugate()).real
    row = {
        "state": tag,
        "frequency_Hz": frequency,
        "phase_sequence": sequence,
        "current_amplitude_A": amplitude,
        "conductivity_MSm_m": conductivity,
        "mesh_mm": mesh,
        "boundary_radius_mm": radius,
        "elements": mo_numelements(),
        "Fx_N": fx,
        "Fy_N": fy,
        "lorentz_x_N": lorentz,
        "loss_W": loss,
        "summed_plate_current_real_A": plate_current.real,
        "summed_plate_current_imag_A": plate_current.imag,
        "summed_plate_current_abs_A": abs(plate_current),
        "circuit_active_power_W": active,
        "model": tag + ".fem",
    }
    ROWS.append(row)
    current_residual = abs(plate_current) / max(amplitude, 1e-15)
    check(
        "Изолированная пластина: нулевой суммарный ток по интегралу J: " + tag,
        current_residual,
        1e-6,
        current_residual < 1e-6,
    )
    if conductivity:
        check(
            "Баланс активной мощности: " + tag,
            relative(active, loss),
            0.005,
            relative(active, loss) < 0.005,
        )
        check(
            "Тяга: напряжения Максвелла и сила Лоренца: " + tag,
            relative(fx, lorentz),
            0.025,
            relative(fx, lorentz) < 0.025,
        )
    mo_close()
    mi_close()
    return row


mesh_rows = [
    field("сетка-" + str(k + 1), mesh) for k, mesh in enumerate((1.2, 0.8, 0.5))
]
fine = mesh_rows[-1]
check(
    "Сходимость тяги на двух последних сетках",
    relative(mesh_rows[-2]["Fx_N"], fine["Fx_N"]),
    0.025,
    relative(mesh_rows[-2]["Fx_N"], fine["Fx_N"]) < 0.025,
)
check(
    "Сходимость потерь на двух последних сетках",
    relative(mesh_rows[-2]["loss_W"], fine["loss_W"]),
    0.025,
    relative(mesh_rows[-2]["loss_W"], fine["loss_W"]) < 0.025,
)
reverse = field("обратное-чередование", 0.5, sequence=-1)
check(
    "Реверс фаз меняет знак тяги",
    relative(reverse["Fx_N"], -fine["Fx_N"]),
    0.01,
    relative(reverse["Fx_N"], -fine["Fx_N"]) < 0.01 and abs(fine["Fx_N"]) > 1e-5,
)
check(
    "Реверс сохраняет джоулевы потери",
    relative(reverse["loss_W"], fine["loss_W"]),
    0.01,
    relative(reverse["loss_W"], fine["loss_W"]) < 0.01,
)
doubled = field("двойной-ток", 0.5, amplitude=200)
check(
    "Тяга пропорциональна квадрату тока",
    relative(doubled["Fx_N"], 4 * fine["Fx_N"]),
    0.001,
    relative(doubled["Fx_N"], 4 * fine["Fx_N"]) < 0.001,
)
check(
    "Потери пропорциональны квадрату тока",
    relative(doubled["loss_W"], 4 * fine["loss_W"]),
    0.001,
    relative(doubled["loss_W"], 4 * fine["loss_W"]) < 0.001,
)
insulator = field("непроводящая-пластина", 0.5, conductivity=0)
check(
    "Без проводимости интеграл силы Лоренца равен нулю",
    abs(insulator["lorentz_x_N"]),
    1e-12,
    abs(insulator["lorentz_x_N"]) < 1e-12,
)
check(
    "Без проводимости джоулевы потери равны нулю",
    abs(insulator["loss_W"]),
    1e-12,
    abs(insulator["loss_W"]) < 1e-12,
)
zero_residual = abs(insulator["Fx_N"]) / max(abs(fine["Fx_N"]), 1e-15)
check(
    "Остаток интеграла напряжений в воздухе относительно рабочей тяги",
    zero_residual,
    1e-4,
    zero_residual < 1e-4,
)
far = field("удалённая-граница", 0.5, radius=600)
check(
    "Влияние удаления внешней границы на тягу",
    relative(far["Fx_N"], fine["Fx_N"]),
    0.02,
    relative(far["Fx_N"], fine["Fx_N"]) < 0.02,
)
for frequency in (10, 25, 100, 200, 500):
    field("частота-" + str(frequency), 0.5, frequency=frequency)
finish({"x": "frequency_Hz", "y": "Fx_N"})