"""Законченный уменьшенный ENSEM-подобный опыт: ток → потери → температура.
FEMM 5/MicroPython. Диск r=20 мм, t=2 мм; mu_r=100, sigma=1.67 МСм/м.
Катушка r=5..18 мм, z=-5..-3 мм, N=10, I_peak=sqrt(2)*10 А; тонкие витки.
f=1 и 25 кГц; сетки delta/3 и delta/6, затем Small Skin Depth на том же контуре.
Тепловая модель: k=15 Вт/(м К), h=10 Вт/(м² К), Ta=293.15 К; потери при 25 кГц.
Потери усредняются по объёму с сохранением мощности; стационарная температура
не предсказывает время закипания и не учитывает воду/движение/зависимость sigma(T).
"""

# Каталог сценария определяется явно: запуск не зависит от чужого 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-induction-heating")
rows = []
for frequency, mode, divisor in [
    (1000, "Bulk", 3),
    (1000, "Bulk", 6),
    (1000, "Impedance", 3),
    (25000, "Bulk", 3),
    (25000, "Bulk", 6),
    (25000, "Impedance", 3),
]:
    delta = math.sqrt(2 / (2 * math.pi * frequency * MU0 * 100 * 1.67e6))
    newdocument(0)
    mi_probdef(frequency, "millimeters", "axi", 1e-9, 1, 30)
    material("Air")
    material("Metal", mu=100, sigma=1.67)
    mi_addcircprop("Coil", math.sqrt(2) * 10, 1)
    mi_drawrectangle(5, -5, 18, -3)
    label(10, -4, "Air", 0.4, "Coil", 10, 1)
    mi_drawrectangle(0, 0, 20, 2)
    label(25, 0, "Air", 1.5)
    mi_makeABC(7, 70, 0, 0, 0)
    if mode == "Bulk":
        label(10, 1, "Metal", delta * 1000 / divisor, group=2)
    else:
        mi_selectsegment(0, 1)
        mi_deleteselectedsegments()
        mi_addboundprop("SSD", 0, 0, 0, 0, 100, 1.67, 0, 0, 1)
        for x, y in [(10, 0), (20, 1), (10, 2)]:
            mi_selectsegment(x, y)
        mi_setsegmentprop("SSD", delta * 1000 / 3, 0, 0, 2)
        mi_clearselected()
    solve(folder, "%dHz-%s-%d" % (frequency, mode, divisor))
    current, voltage, linkage = mo_getcircuitproperties("Coil")
    transmitted = -(2 * math.pi * frequency * linkage * current.conjugate()).imag / 2
    row = {
        "frequency_Hz": frequency,
        "mode": mode,
        "mesh_divisor": divisor,
        "delta_mm": delta * 1000,
        "P_from_linkage_W": transmitted,
        "elements": mo_numelements(),
    }
    if mode == "Bulk":
        mo_groupselectblock(2)
        row["P_volume_W"] = mo_blockintegral(4).real
    rows.append(row)
    close_magnetic()
    save_json(folder, {"magnetic": rows})
power = rows[4]["P_volume_W"]
volume = math.pi * 0.02**2 * 0.002
newdocument(2)
hi_probdef("millimeters", "axi", 1e-9, 1, 30)
hi_addmaterial("Metal", 15, 15, power / volume, 3.9)
hi_addboundprop("Convection", 2, 0, 0, 293.15, 10, 0)
hi_drawrectangle(0, 0, 20, 2)
hi_addblocklabel(10, 1)
hi_selectlabel(10, 1)
hi_setblockprop("Metal", 0, 0.25, 1)
hi_clearselected()
for x, y in [(10, 0), (20, 1), (10, 2)]:
    hi_selectsegment(x, y)
hi_setsegmentprop("Convection", 0.25, 0, 0, 0)
hi_clearselected()
hi_saveas(os.path.join(folder, "temperature.feh"))
hi_analyze(1)
hi_loadsolution()
temperature = ho_getpointvalues(0, 1)[0]
ho_close()
hi_close()
surface = 2 * math.pi * 0.02**2 + 2 * math.pi * 0.02 * 0.002
data = {
    "magnetic": rows,
    "thermal": {
        "input_W": power,
        "volume_m3": volume,
        "uniform_q_W_m3": power / volume,
        "center_T_K": temperature,
        "lumped_T_K": 293.15 + power / (10 * surface),
    },
}
save_json(folder, data)
assert power > 0 and temperature > 293.15
assert abs(power / rows[4]["P_from_linkage_W"] - 1) < 0.03
print(data)