"""Два длинных круглых провода: Лоренц, WST и контур Максвелла.
Радиус 1 мм, оси x=+-10 мм, I=100 А, глубина 100 мм, DC.
Одинаковые токи притягиваются; сила левого провода направлена по +x.
Концевые эффекты не учитываются. Сетка и внешний радиус проверяются раздельно.
"""

# Каталог сценария определяется явно: запуск не зависит от чужого 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-conductor-force")
rows = []
exact = MU0 * 100**2 * 0.1 / (2 * math.pi * 0.02)
for mesh, radius in [(0.3, 100), (0.15, 100), (0.15, 200)]:
    newdocument(0)
    mi_probdef(0, "millimeters", "planar", 1e-10, 100, 30)
    material("Air")
    mi_addcircprop("I", 100, 1)
    for x, group in [(-10, 1), (10, 2)]:
        circle(x, 0, 1)
        label(x, 0, "Air", mesh, "I", 1, group)
    label(0, 0, "Air", 2 * mesh)
    # Вдали от проводов достаточно редкой сетки; уточняем рабочую область.
    circle(0, 0, 25)
    label(50, 0, "Air", 5)
    a0_circle(radius)
    solve(folder, "mesh-%g-R-%g" % (mesh, radius))
    mo_groupselectblock(1)
    lorentz = [mo_blockintegral(11), mo_blockintegral(12)]
    weighted = [mo_blockintegral(18), mo_blockintegral(19)]
    mo_clearblock()
    contours = []
    for r in [2, 3, 4]:
        mo_clearcontour()
        # Обход по часовой стрелке соответствует внешней нормали FEMM.
        for x, y in [
            (-10 - r, -r),
            (-10 - r, r),
            (-10 + r, r),
            (-10 + r, -r),
            (-10 - r, -r),
        ]:
            mo_addcontour(x, y)
        contours.append([r, list(mo_lineintegral(3))])
    rows.append(
        {
            "mesh_mm": mesh,
            "outer_radius_mm": radius,
            "exact_Fx_N": exact,
            "Lorentz_N": lorentz,
            "WST_N": weighted,
            "Maxwell_contours_N": contours,
            "elements": mo_numelements(),
        }
    )
    close_magnetic()
    save_json(folder, rows)
assert abs(rows[-1]["Lorentz_N"][0] / exact - 1) < 0.02
print(rows)