"""Дифференциальная индуктивность из независимых магнитостатических решений.

Запустите во встроенной Python-консоли FEMM 5, изменив параметры ниже.
Внешние Python, pyFEMM, TCP-подключение и установка пакетов не требуются.
Укажите абсолютные пути с разделителями / и новое имя выходного каталога.
MODELS содержит имена входных моделей для двух уровней сетки. В уточнённой
модели уменьшите положительные размеры элементов у меток вдвое; геометрию,
материалы, витки, источники и граничные условия сохраните. С одной сеткой
расчёт возможен, но сравнение пространственной дискретизации отсутствует.

Вход: VARY — имя изменяемой цепи; READ_CIRCUITS — имена считываемых цепей;
I0 — ток рабочей точки, А; DELTA — положительное приращение тока, А.
Остальные токи берутся из каждой исходной модели и при вариации сохраняются.
Модели должны быть магнитостатическими, с Frequency = 0 и пустым Previous
Solution: требуются полные нелинейные решения, а не Incremental или Frozen.

Выход: samples.json содержит токи, А, потокосцепления, Вб·виток, энергии, Дж,
и размеры сетки. Исторический ключ linkage_Wb обозначает потокосцепление Ψ,
а не магнитный поток отдельного витка. В summary.json записаны производные
по току, Гн, для двух шагов и двух сеток. Завершение расчёта не устанавливает
точность: отдельно сопоставьте чувствительность к шагу и к сгущению сетки.

Ориентация матрицы: строка i — цепь из READ_CIRCUITS, столбец j — VARY;
результат равен L_ij = ∂Ψ_i/∂I_j при фиксированных остальных токах. Этот запуск
даёт один столбец. Для полной матрицы повторите запуск с другим VARY и новым
выходным каталогом, сохраняя один и тот же полный вектор токов рабочей точки.
Знаки взаимных членов зависят от направлений токов и витков входных моделей.
"""

import json
import os

# 1. Настройка: исходные модели открываются заново для каждой точки.
# Пути абсолютные; OUTPUT должен быть новым каталогом с существующим родителем.
MODELS = {"base": "C:/FEMM-lab/nonlinear-3.fem",
          "refined": "C:/FEMM-lab/nonlinear-3-fine.fem"}
OUTPUT = "C:/FEMM-lab/inductance-new-run"
VARY = "Primary"
READ_CIRCUITS = ["Primary", "Secondary"]
I0 = 3.0
DELTA = 0.1

assert DELTA > 0, "Шаг тока должен быть положительным"
assert VARY in READ_CIRCUITS, "Считывайте также возбуждаемую цепь"

# Если каталог существует, mkdir завершится ошибкой: результаты сохраняются.
os.mkdir(OUTPUT)
rows = []
summary = {}

# 2. Для каждой сетки рассчитываем пять токов: I0 ± Δ, I0 ± Δ/2 и I0.
# Общие крайние и центральные точки позволяют сравнить два шага производной
# без повторных решений одинакового режима. Порядок samples задаётся offsets.
for mesh in MODELS:
    samples = []
    offsets = [-DELTA, -DELTA/2, 0.0, DELTA/2, DELTA]
    for number in range(len(offsets)):
        # До открытия проверяем отсутствие ссылки на предыдущее решение.
        # Нужны полные независимые решения, а не Incremental или Frozen.
        with open(MODELS[mesh], "r") as stream:
            for line in stream.read().split("\n"):
                if line.startswith("[PrevSoln]"):
                    assert line.split("=", 1)[1].strip() == '""', "Очистите Previous Solution в исходной модели"

        # Команды mi относятся к активному входному магнитному документу.
        # info[1] — частота, Гц; координаты и глубина наследуются из модели.
        opendocument(MODELS[mesh])
        info = mi_getprobleminfo()
        assert abs(info[1]) < 1e-12, "Нужна магнитостатическая модель (Frequency = 0)"

        # Сначала сохраняем рабочую копию, затем меняем ток цепи VARY.
        # Индекс свойства 1 означает ток, А. Остальные цепи не изменяются.
        path = OUTPUT + "/" + mesh + "-" + str(number) + ".fem"
        mi_saveas(path)
        mi_modifycircprop(VARY, 1, I0 + offsets[number])
        mi_saveas(path)
        mi_analyze()
        mi_loadsolution()

        # mo читает только загруженное решение этой рабочей копии.
        # В статике потокосцепление должно быть действительным; мнимая часть
        # проверяется до преобразования, чтобы не скрыть неверную постановку.
        row = {"mesh": mesh, "offset_A": offsets[number], "circuits": {}}
        for name in READ_CIRCUITS:
            current, voltage, linkage = mo_getcircuitproperties(name)
            assert abs(complex(linkage).imag) < 1e-10, "Ожидалось действительное DC-потокосцепление"
            row["circuits"][name] = {"current_A": complex(current).real,
                                      "linkage_Wb": complex(linkage).real}

        # Снимаем предыдущее выделение и выбираем все блоки для полной энергии.
        # Интеграл 2 — энергия W, интеграл 17 — коэнергия W', оба в джоулях.
        mo_clearblock()
        mo_groupselectblock()
        row["energy_J"] = complex(mo_blockintegral(2)).real
        row["coenergy_J"] = complex(mo_blockintegral(17)).real
        row["nodes"] = int(mo_numnodes())
        row["elements"] = int(mo_numelements())
        samples.append(row)
        rows.append(row)

        # Закрываем результат и вход перед независимым расчётом следующего тока.
        mo_close()
        mi_close()

    # 3. Для каждой считываемой цепи строим одну строку столбца L[:, VARY].
    # Центральные разности используют симметричные точки, что устраняет
    # первый порядок ошибки по шагу для гладкой зависимости Ψ_i(I_j).
    entry = {}
    for name in READ_CIRCUITS:
        psi = [sample["circuits"][name]["linkage_Wb"] for sample in samples]

        # Полный шаг Δ: знаменатель 2Δ; половинный шаг Δ/2: знаменатель Δ.
        coarse_step = (psi[4]-psi[0])/(2*DELTA)
        half_step = (psi[3]-psi[1])/DELTA
        entry[name] = {"Psi_at_bias_Wb": psi[2], "Ldynamic_delta_H": coarse_step,
                       "Ldynamic_halfdelta_H": half_step,
                       "step_change_relative": abs(half_step-coarse_step)/max(abs(half_step), 1e-30)}

    # Отношения Ψ/I и 2W/I² сохранены как дополнительные диагностические числа.
    # В нелинейной или многоконтурной системе они в общем случае не равны
    # дифференциальной самоиндуктивности. Деление на нулевой I0 исключается.
    if I0 != 0:
        entry[VARY]["Psi_over_I_H"] = samples[2]["circuits"][VARY]["linkage_Wb"]/I0
        entry[VARY]["twoW_over_I2_H"] = 2*samples[2]["energy_J"]/(I0*I0)
    entry["elements"] = samples[2]["elements"]
    summary[mesh] = entry

# 4. Записываем первичные отсчёты отдельно от производных: это позволяет
# повторить обработку без нового полевого расчёта и проследить каждый результат.
with open(OUTPUT + "/samples.json", "w") as stream:
    json.dump(rows, stream)
with open(OUTPUT + "/summary.json", "w") as stream:
    json.dump(summary, stream)
print(json.dumps(summary))
print("Применяйте отношения Ψ/I и 2W/I² только при предпосылках, указанных в уроке.")
print("Сравните оба шага и обе сетки: завершённый расчёт ещё не подтверждает точность.")