"""Точный эталон Hall и поисковой катушки в однородном поле FEMM 5.

Среда: встроенная Python-консоль FEMM 5; внешние пакеты не нужны.
Вход: B0 в Т, размеры и шаг сетки в мм, глубина Planar в мм, число витков.
Выход: четыре новые модели .fem и results.json с B [Т], H [А/м], Phi [Вб].
OUTPUT должен отсутствовать: при повторе выберите новое имя каталога.
Поле задаётся граничным потенциалом в воздухе; источник-магнит не моделируется.
Hall и катушка — идеальные невозмущающие наблюдатели без корпуса и калибровки.
"""

import json
import math
import os

OUTPUT = "measurement-field-run"
# mkdir намеренно прерывает повторный запуск, сохраняя предыдущие результаты.
os.mkdir(OUTPUT)

MU0 = 4 * math.pi * 1e-7  # Учебное приближение, Гн/м.
B0 = 0.01
WIDTH_MM = 20.0
HEIGHT_MM = 20.0
DEPTH_MM = 10.0
AREA = WIDTH_MM * DEPTH_MM * 1e-6  # мм² -> м²; поверхность потока x-z.
TURNS = 100
FREQ_HZ = 50.0


def real(value):
    """Извлечь действительную часть результата магнитостатического расчёта."""
    z = complex(value)
    # Комплексная компонента здесь указывает на неподходящий режим расчёта.
    assert abs(z.imag) < 1e-10, "Нужна магнитостатика"
    return z.real


# Две сетки и два добавочных постоянных потенциала проверяют сходимость
# и калибровочную инвариантность: прибавка к A_z не должна менять B.
rows = []
for mesh_mm in (0.25, 0.125):
    for a_offset in (0.0, 0.0005):
        newdocument(0)
        mi_probdef(0, "millimeters", "planar", 1e-10, DEPTH_MM, 25)
        mi_addmaterial("Air-reference", 1, 1, 0, 0, 0)

        # В Planar: B_y = -dA_z/dx. A_z измеряется в Вб/м,
        # поэтому ширину для разности потенциалов переводим из мм в м.
        mi_addboundprop("Left-A", a_offset, 0, 0, 0, 0, 0, 0, 0, 0)
        right_a = a_offset - B0 * WIDTH_MM * 1e-3
        mi_addboundprop("Right-A", right_a, 0, 0, 0, 0, 0, 0, 0, 0)
        corners = ((0, 0), (WIDTH_MM, 0), (WIDTH_MM, HEIGHT_MM), (0, HEIGHT_MM))
        for x, y in corners:
            mi_addnode(x, y)
        for k in range(4):
            x1, y1 = corners[k]
            x2, y2 = corners[(k + 1) % 4]
            mi_addsegment(x1, y1, x2, y2)
        for x, name in ((0, "Left-A"), (WIDTH_MM, "Right-A")):
            mi_selectsegment(x, HEIGHT_MM / 2)
            mi_setsegmentprop(name, 0, 1, 0, 0)
            mi_clearselected()

        # Верх и низ имеют естественное условие dA_z/dn = 0.
        # Отключённый автоматический размер задаёт независимый шаг сетки.
        mi_addblocklabel(WIDTH_MM / 2, HEIGHT_MM / 2)
        mi_selectlabel(WIDTH_MM / 2, HEIGHT_MM / 2)
        mi_setblockprop("Air-reference", 0, mesh_mm, "<None>", 0, 1, 0)
        mi_clearselected()
        stem = "mesh-" + str(mesh_mm) + "-gauge-" + str(a_offset)
        mi_saveas(OUTPUT + "/" + stem + ".fem")
        mi_analyze(1)
        mi_loadsolution()
        mo_smooth("off")

        # Позиции 1, 2, 5, 6 в mo_getpointvalues — Bx, By, Hx, Hy.
        point = mo_getpointvalues(WIDTH_MM / 2, HEIGHT_MM / 2)
        bx, by, hx, hy = [real(point[k]) for k in (1, 2, 5, 6)]

        # Потоковый интеграл уже учитывает глубину Planar: повторный
        # множитель DEPTH_MM недопустим. Знак задаёт направление контура.
        mo_clearcontour()
        mo_addcontour(0, HEIGHT_MM / 2)
        mo_addcontour(WIDTH_MM, HEIGHT_MM / 2)
        integral = mo_lineintegral(0)
        signed_flux = real(integral[0])
        average_bn = real(integral[1])
        mo_clearcontour()
        mo_addcontour(WIDTH_MM, HEIGHT_MM / 2)
        mo_addcontour(0, HEIGHT_MM / 2)
        reversed_flux = real(mo_lineintegral(0)[0])

        # Интеграл H·dl вдоль положительного y на пути длиной 16 мм, А.
        mo_clearcontour()
        mo_addcontour(WIDTH_MM / 2, 2)
        mo_addcontour(WIDTH_MM / 2, 18)
        magnetic_voltage = real(mo_lineintegral(1)[0])
        row = {
            "mesh_mm": mesh_mm,
            "a_offset_Wb_per_m": a_offset,
            "nodes": int(mo_numnodes()),
            "elements": int(mo_numelements()),
            "Bx_T": bx,
            "By_T": by,
            "Hx_A_per_m": hx,
            "Hy_A_per_m": hy,
            "signed_flux_Wb": signed_flux,
            "average_Bn_T": average_bn,
            "reversed_flux_Wb": reversed_flux,
            "magnetic_voltage_A": magnetic_voltage,
            # Идеальный Hall: нулевой уровень 2,5 В и чувствительность 20 В/Т.
            "Hall_output_V": 2.5 + 20 * by,
            "linkage_magnitude_Wb": TURNS * abs(signed_flux),
            # DC-решение задаёт амплитуду предполагаемой синусоиды потока.
            # Это аналитическая оценка ЭДС, а не временное решение FEMM.
            "emf_rms_50Hz_V": (
                2 * math.pi * FREQ_HZ * TURNS * abs(signed_flux) / math.sqrt(2)
            ),
        }

        # Допуски заданы в единицах соответствующих величин и сопоставлены
        # с точным решением: B_y = B0, H_y = B0/mu0, |Phi| = B0*S.
        assert abs(bx) < 1e-8 and abs(by - B0) < 1e-8
        assert abs(hx) < 0.01 and abs(hy - B0 / MU0) < 0.01
        assert abs(abs(signed_flux) - B0 * AREA) < 1e-10
        assert abs(abs(average_bn) - B0) < 1e-8
        assert abs(signed_flux + reversed_flux) < 1e-12
        assert abs(magnetic_voltage - B0 / MU0 * 0.016) < 0.001
        rows.append(row)
        mo_close()
        mi_close()

# Каждая пара строк относится к одной сетке и двум значениям прибавки к A_z.
for k in (0, 2):
    assert abs(rows[k]["By_T"] - rows[k + 1]["By_T"]) < 1e-8
    assert abs(rows[k]["signed_flux_Wb"] - rows[k + 1]["signed_flux_Wb"]) < 1e-10
assert rows[2]["elements"] > rows[0]["elements"], "Нужны разные сетки"

# Сохраняем численные величины вместе с происхождением и границами модели.
result = {
    "origin": "native_FEMM_static_uniform_air_field",
    "scope": "Point B/H, oriented flux, 16 mm H path, mesh and gauge checks; ideal passive sensors; no physical calibration or device validation.",
    "input": {"B_T": B0, "width_mm": WIDTH_MM, "depth_mm": DEPTH_MM, "turns": TURNS},
    "expected": {"flux_magnitude_Wb": B0 * AREA, "Hy_A_per_m": B0 / MU0},
    "rows": rows,
    "passed": len(rows) == 4,
}
with open(OUTPUT + "/results.json", "w") as stream:
    json.dump(result, stream)
print(json.dumps(result))