"""Независимый тест транспортного скин-эффекта FEMM5, не пересказ вывода2020г.
Круглая медь a=1мм, длина100мм, sigma=58МСм/м,I_peak=100А. f=1 и100кГц.
Точный внутренний импеданс через ряды J0,J1; внешняя индуктивность исключается.
Две сетки на каждой частоте; реальные провода, LamType=0, задан именно ток цепи.
"""

# Каталог сценария определяется явно: запуск не зависит от чужого CWD.
import os
import sys

os.chdir(os.path.dirname(os.path.abspath(sys.argv[0])))
exec(open("common.py").read())
folder = result_folder("results-skin-benchmark")


def bessel(order, z):
    term = 1 + 0j if order == 0 else z / 2
    total = term
    for k in range(1, 200):
        term *= -(z * z / 4) / (k * (k + order))
        total += term
        if abs(term) < 1e-14 * max(abs(total), 1):
            return total
    raise RuntimeError("Ряд Бесселя не сошёлся")


rows = []
for frequency in [1000, 100000]:
    omega = 2 * math.pi * frequency
    delta = math.sqrt(2 / (omega * MU0 * 58e6))
    k = (1 - 1j) / delta
    exact = (
        0.1
        * k
        * bessel(0, k * 0.001)
        / (2 * math.pi * 0.001 * 58e6 * bessel(1, k * 0.001))
    )
    for divisor in [3, 6]:
        newdocument(0)
        mi_probdef(frequency, "millimeters", "planar", 1e-10, 100, 30)
        material("Air")
        material("Cu", sigma=58)
        mi_addcircprop("Wire", 100, 1)
        circle(0, 0, 1, 1)
        label(0, 0, "Cu", min(0.3 / divisor, delta * 1000 / divisor), "Wire", 1, 1)
        circle(0, 0, 3, 1)
        label(2, 0, "Air", 0.15 / divisor)
        label(5, 0, "Air", 0.3)
        a0_circle(10)
        solve(folder, "skin-%d-%d" % (frequency, divisor))
        current, voltage, flux = mo_getcircuitproperties("Wire")
        mo_groupselectblock(1)
        power = mo_blockintegral(4).real
        internal = voltage / current - 1j * omega * MU0 * 0.1 / (
            2 * math.pi
        ) * math.log(10)
        rows.append(
            {
                "frequency_Hz": frequency,
                "divisor": divisor,
                "delta_mm": delta * 1000,
                "Z_internal_exact_ohm": pair(exact),
                "Z_internal_FE_ohm": pair(internal),
                "R_from_loss_ohm": 2 * power / 100**2,
                "P_FE_W": power,
                "P_exact_W": 100**2 * exact.real / 2,
                "relative_R_error": abs(2 * power / 100**2 / exact.real - 1),
            }
        )
        close_magnetic()
        save_json(folder, rows)
assert max(row["relative_R_error"] for row in rows) < 0.03
print(rows)