"""Проводящая труба: FE и независимое интегрирование точного радиального ОДУ.
FEMM5/MicroPython, без scipy. ri=12.7мм,ro=25.4мм,mu_r=1000,sigma=10МСм/м,
f=1Гц,I_peak=100А, глубина1м. Сначала труба I=0, затем заземлённая на бесконечности.
RK4 решает линейное комплексное ОДУ; 2000 и 4000 шагов проверяют интегрирование.
"""

# Каталог сценария определяется явно: запуск не зависит от чужого 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-tube")
ri, ro = 0.0127, 0.0254
alpha2 = 1j * 2 * math.pi * 1000 * MU0 * 1e7


def integrate(h_initial, d_initial, count):
    step = (ro - ri) / count
    h, d = complex(h_initial), complex(d_initial)
    values = [[ri, h, d]]

    def slope(r, h, d):
        return d, (1 / r**2 + alpha2) * h - d / r

    for n in range(count):
        r = ri + n * step
        a, b = slope(r, h, d)
        c, e = slope(r + step / 2, h + a * step / 2, d + b * step / 2)
        f, g = slope(r + step / 2, h + c * step / 2, d + e * step / 2)
        k, l = slope(r + step, h + f * step, d + g * step)
        h += step * (a + 2 * c + 2 * f + k) / 6
        d += step * (b + 2 * e + 2 * g + l) / 6
        values.append([r + step, h, d])
    return values


def reference(count):
    a = integrate(100 / (2 * math.pi * ri), 0, count)
    b = integrate(0, 1, count)
    coefficient = (100 / (2 * math.pi * ro) - a[-1][1]) / b[-1][1]
    return [
        [v[0], v[1] + coefficient * w[1], v[2] + coefficient * w[2]]
        for v, w in zip(a, b)
    ]


ref = reference(4000)
ref_coarse = reference(2000)
ode_difference = max(abs(ref[2 * n][1] - ref_coarse[n][1]) for n in range(2001))
rows = []
for constrained, mesh in [(True, 0.8), (True, 0.4), (False, 0.4)]:
    newdocument(0)
    mi_probdef(1, "millimeters", "planar", 1e-10, 1000, 30)
    material("Air")
    material("Steel", 1000, 10)
    mi_addcircprop("Wire", 100, 1)
    mi_addcircprop("Tube", 0, 0)
    circle(0, 0, 2.54)
    circle(0, 0, 12.7)
    a0_circle(25.4)
    label(0, 0, "Air", mesh, "Wire", 1, 1)
    label(5, 0, "Air", mesh)
    label(20, 0, "Steel", mesh, "Tube" if constrained else "<None>", 1, 2)
    solve(folder, "Izero-%s-mesh-%g" % (constrained, mesh))
    samples = []
    numerator, denominator = 0, 0
    for index in range(100, 4000, 200):
        r, h, derivative = ref[index]
        exact_j = derivative + h / r
        p = mo_getpointvalues(r * 1000, 0)
        j = (p[7] + p[8]) * 1e6
        samples.append(
            {
                "r_m": r,
                "H_exact": pair(h),
                "H_FE": pair(p[6]),
                "J_exact": pair(exact_j),
                "J_FE": pair(j),
            }
        )
        numerator += abs(j - exact_j) ** 2
        denominator += abs(exact_j) ** 2
    mo_groupselectblock(2)
    rows.append(
        {
            "Izero": constrained,
            "mesh_mm": mesh,
            "tube_current_A": pair(mo_blockintegral(7)),
            "power_W": mo_blockintegral(4).real,
            "relative_RMS_J_error": math.sqrt(numerator / denominator),
            "ode_mesh_difference_A_per_m": ode_difference,
            "samples": samples,
        }
    )
    close_magnetic()
    save_json(folder, rows)
assert ode_difference < 1e-3
assert abs(complex(*rows[1]["tube_current_A"])) < 1e-5
assert rows[1]["relative_RMS_J_error"] < 0.08
print(rows)