"""Нативный эталон коаксиального конденсатора с тремя уровнями сетки.

Запуск: встроенный Python FEMM 5, отдельный пустой текущий рабочий каталог.
Код создаёт модели coax-*.fee, их решения и coax-results.json с фиксированными
именами. Сохраните открытые документы до запуска. Внешние пакеты не нужны.

Вход: радиусы a и b, мм; напряжение voltage, В; относительная проницаемость er;
осевая длина 100 мм задаётся глубиной плоской задачи. Это поперечное сечение
длинного коаксиального конденсатора без торцевого поля, а не осесимметричная
постановка FEMM. Диэлектрик линейный, однородный, потерь и свободного заряда нет.

Алгоритм: для каждой пары размера элемента и углового шага заново построить
модель, решить её и сравнить энергию, потенциал и радиальное поле с формулами.
Выход: ёмкость, Ф; потенциал, В; напряжённость, В/м; относительные ошибки.
Сгущение одновременно меняет сетку и аппроксимацию дуг; критерии ниже относятся
к этому эталону и не устанавливают общий порядок сходимости решателя.
"""

import math, json

out = {"geometry": "planar coaxial cross-section; axial length 100 mm", "cases": [], "passed": False}

# 1. Исходные данные. В ln(b/a) входят одинаковые единицы радиусов;
# множитель 0,1 в формуле ёмкости — осевая длина в метрах.
a = 5.0
b = 15.0
voltage = 1000.0
er = 3.2
eps0 = 8.8541878128e-12
Cexact = 2 * math.pi * eps0 * er * .1 / math.log(b / a)

# Геометрическое среднее радиусов даёт V(r) = voltage/2 при V(b) = 0.
# Точка (r, 0) лежит на положительной оси x: E_x совпадает с E_r > 0.
r = (a * b)**.5

# 2. Три независимых построения: mesh в мм, angle в градусах.
for mesh, angle in [(0.8, 2.0), (0.4, 1.0), (0.2, 0.5)]:
    newdocument(1)
    ei_probdef("millimeters", "planar", 1e-10, 100, 30)
    ei_addmaterial("Insulation", er, er, 0)
    ei_addboundprop("Inner", voltage, 0, 0, 0, 0)
    ei_addboundprop("Outer", 0, 0, 0, 0, 0)

    # Каждая окружность состоит из двух полуокружностей; выбранным дугам
    # назначаются потенциал электрода и максимальный угловой шаг сегмента.
    for radius, name in [(a, "Inner"), (b, "Outer")]:
        ei_addnode(radius, 0)
        ei_addnode(-radius, 0)
        ei_addarc(radius, 0, -radius, 0, 180, angle)
        ei_addarc(-radius, 0, radius, 0, 180, angle)
        ei_selectarcsegment(0, radius)
        ei_selectarcsegment(0, -radius)
        ei_setarcsegmentprop(angle, name, 0, 0, "<None>")
        ei_clearselected()

    # Диэлектрик занимает кольцо. Внутренний электрод исключается из сетки;
    # область за внешним электродом в эту замкнутую постановку не включена.
    ei_addblocklabel(r, 0)
    ei_selectlabel(r, 0)
    ei_setblockprop("Insulation", 0, mesh, 0)
    ei_clearselected()
    ei_addblocklabel(0, 0)
    ei_selectlabel(0, 0)
    ei_setblockprop("<No Mesh>", 1, 0, 0)
    ei_clearselected()

    # Сохранение предшествует решению; команды eo читают только загруженный
    # результат текущей модели. У каждого уровня сетки своё имя файла.
    filename = "coax-" + str(mesh) + ".fee"
    ei_saveas(filename)
    ei_analyze(1)
    ei_loadsolution()
    values = eo_getpointvalues(r, 0)

    # Выбираем весь диэлектрик: интеграл 0 возвращает энергию W, Дж.
    # Индекс [0] извлекает значение из возвращаемой встроенным API коллекции.
    eo_selectblock(r, 0)
    energy = eo_blockintegral(0)[0]
    capacitance = 2 * energy / voltage**2

    # Радиус переводится из мм в м для E_r = U/[r ln(b/a)], В/м.
    # В eo_getpointvalues: [0] — V, В; [3] — E_x, В/м.
    e_exact = voltage / ((r * .001) * math.log(b / a))
    out["cases"].append({
        "mesh_mm": mesh,
        "arc_step_deg": angle,
        "model": filename,
        "capacitance_F": capacitance,
        "exact_C_F": Cexact,
        "C_relative_error": abs(capacitance - Cexact) / Cexact,
        "V_mid": values[0],
        "E_mid": values[3],
        "exact_E_mid": e_exact,
        "E_relative_error": abs(values[3] - e_exact) / e_exact,
        "V_relative_error": abs(values[0] - voltage / 2) / (voltage / 2),
    })
    eo_close()
    ei_close()

# 3. Проверяем последнюю сетку: ошибка C < 0,2 %, E < 1 %, V < 0,2 %;
# дополнительно ошибка ёмкости должна снизиться относительно первой сетки.
# Промежуточная сетка сохраняется в журнале для самостоятельного сравнения.
out["passed"] = (
    out["cases"][-1]["C_relative_error"] < .002
    and out["cases"][-1]["E_relative_error"] < .01
    and out["cases"][-1]["V_relative_error"] < .002
    and out["cases"][-1]["C_relative_error"] < out["cases"][0]["C_relative_error"]
)
# Сохраняем первичные значения каждой сетки вместе с аналитическими эталонами;
# флаг passed относится к перечисленным проверкам, а не к любой геометрии.
with open("coax-results.json", "w") as f:
    json.dump(out, f)
print(json.dumps(out))
assert out["passed"]